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This article reviews lattice QCD results for the light hadron spectrum. We give an 
overview of different formulations of lattice QCD, with discussions on the fermion dou- 
bling problem and improvement programs. We summarize recent developments in al- 
gorithms and analysis techniques, that render calculations with light, dynamical quarks 
feasible on present day computer resources. Finally, we summarize spectrum results for 
ground state hadrons and resonances using various actions. 
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I. INTRODUCTION 



Already at the beginning of the 19 century, it was 



speculated by (Prout 1815) that the hydrogen atom was 



the basic building block for all other atoms. The mass 
of the proton, as mass of the hydrogen atom, was known 
within a factor 2 accuracy already during the late 19"^ 



century (Loschmidt 



1865). Later, the development of 



mass spectrometry (Goldstein 1886) allowed a precision 



measurement of the e/m ratio of the hydrogen nucleus 
(Thomson 1907 Wien 1902 1 and following the discov- 



ery of the atomic nucleus by (Rutherford 1911 ), he could 



show that hydrogen nuclei were present in other nuclei 



(Rutherford 1919) and coined for them the name pro- 



tons. 



The neutron was discovered 13 years later by (Chad 



wick 1932), who also determined its mass with a 2 per 



mil accuracy. The first meson to be discovered was the 



pion (Lattes et al. 1947), shortly followed by the kaon 



(Rochester and Butler 1947) and the A (Seriff et al. 



1950), the first strange particles. While these discoveries 



were made in cosmic ray experiments, the first resonance. 



the A was discovered by (Brueckner 1952 ) at a cyclotron 



source. During the following years, these modern accel- 
erators lead to a proliferation of hadronic states and it 
became obvious that they could not all be regarded as 
elementary. 

This large number of hadronic states could first be 



successfully described by their quark substructure (Gel! 
Mann [1961 ), for which finally quantum chromodynamics 
(QCD) was found as the dynamical theory by (Fritzsch 



et al. 1973). With the discovery of asymptotic free- 



dom (Gross and Wilczek 1973 Politzer 1973), which 



built on earlier work regarding the renormalizability of 



IJTJ nonabelian gauge theories by ('t Hooft and Veltman 



1972), and the qualitative understanding of the con- 



finement phenomenon (Wilson 1974) a coherent pic- 
ture of the strong interaction finally emerged. At en- 
ergies that are large compared to the typical QCD scale 
Aqcd ~ 250 MeV the coupling is small and quarks and 
massless gluons emerge as the fundamental degrees of 
freedom. At low energies however, the spectrum of QCD 



consists of quark-gluon bound states that one would like 
to identify with the experimentally observed hadrons. Al- 
though this qualitative picture is quite compelling, it is 
nevertheless very difficult to solve QCD in the low en- 
ergy regime, where it is a strongly coupled theory, and 
predict respectively postdict the hadron spectrum from 
first principles. 

In the present review, we summarize the current state 
of the art of computing the light hadron spectrum, i.e. 
the spectrum of hadrons with exclusively up, down and 
strange valence quarks, directly in QCD. In sect.|n]we in- 
troduce the primary tool to study QCD in the nonpertur- 
bative regime - lattice QCD. We review various possible 
discretizations of continuum QCD in view of their use- 
fulness for ab initio calculations of light hadron masses 
with a small and controlled total uncertainty. In sect. |III[ 
we review current methods of extracting hadron masses 
in lattice QCD. We discuss efficient ways of extracting 
ground state masses as well as current methods to over- 
come the challenges in singlet and excited state spec- 
troscopy. In sect. |IV| we review methods for obtaining 
predictions at the physical point (the point in parame- 
ter space at which the quark masses have their physical 
values) in the infinite volume continuum theory. Finally 
in sect. |V] we summarize present and notable past re- 
sults and conclude with an overview of the present un- 
derstanding of the hadron spectrum from lattice QCD. 
As a convention when quoting lattice results the first er- 
ror is statistical and the second one (if given) systematic 
unless explicitly noted otherwise. 



II. LATTICE TECHNIQUES 

Lattice field theory is in most cases the only known 
systematic way of nonperturbatively computing Greens 
functions in quantum field theories. It is especially use- 
ful in contexts where perturbative treatment is usually 
inadequate, which is the case in low energy QCD. 



Lattice gauge theory was introduced by ( Wilson 



1974 



tar 



and recent overviews include (DeGrand and De- 
2000; Gattringcr and Lang 



2006 Di Pierro 



Gupta[ |1997| [Montvay and Munster, ,1994^ |Rothe| 



2010 



2005 



Smit||2002[ l^ 

In general, a nonperturbative lattice calculation pro- 
ceeds in 3 steps. First, one introduces a UV regulator 
into the theory by means of a finite spacetime lattice. 
Then one computes Greens functions in this discretized 
theory by means of stochastic integration of the path in- 
tegral and finally one removes the regulator in order to 
obtain the continuum result. The last step is possible in 



theories where the coupling does not diverge in the UV 
regime. Due to asymptotic freedom, QCD does belong 
to this class of theories and the cutoff can be removed. 

In this section we will mainly focus on the first step 
of the above procedure, the regularization of QCD on a 
spacetime lattice. This regularization is not unique and 
the ambiguity is reflected in the wide variety of lattice 
regularizations of QCD that are in use today each car- 
rying various advantages and disadvantages. We start 
with a brief introduction to the path integral formalism 
in sect. III. Al and the basics of the lattice discretization of 
QCD in sect. |II.B[ We then introduce the basic concepts 
of the stochastic evaluation of the discretized path inte- 
gral in sect. |II.C| which are necessary to understand the 
further developments of sect. [Tl.D| and sect. |II.E] where we 
discuss how to obtain efficient lattice regularized theories 
that actually go over into QCD in upon removal of the 
cutoff. Finally, in sect. |II.F] we briefly discuss anisotropic 
lattice regularizations of QCD that are relevant for ex- 
cited state spectroscopy. 

The need for an efficient regularization arises particu- 
larly due to the smallness of the light quark masses com- 
pared to the intrinsic QCD scale Aqcq. On the one hand, 
the physical size of the lattice needs to be much larger 
than the correlation length of the system which in turn 
is given by the inverse of the mass of the lightest particle 
in the spectrum, the pion. On the other hand, the lat- 
tice cutoff needs to be much larger than Aqcd in order 
to not miss a substantial fraction of the nonperturbative 
dynamics. These two requirements combined necessitate 
a large number of lattice points if one would like to per- 
form nonperturbative lattice QCD calculations at physi- 
cally light quark masses. In connection with the fermion 
doubling problem the smallness of the light quark masses 
causes yet further problems that are discussed in detail 
in sect. HLD] 

Because of these effects, lattice QCD calculations un- 
til very recently were restricted to quark masses larger 
(and in most cases substantially so) than the physical 
ones. This in turn necessitated an extrapolation in the 
light quark mass to the physical point in addition to 
the already necessary continuum extrapolation. The fact 
that physically light quark masses have been reached 



by reweighting (Aoki et al. 2010) or directly, at large 



volumes and several different values of the cutoff (Durr 
et al. 2011a|b ) is to a large extent due to recent advances 
in the construction of efficient lattice regularizations that 
will be reviewed in this section. 



A. Basics of the Path Integral formalism 



^ Independent developments of Smit and Polyakov were never 
published, see e.g. (Wilson 2005) 



See also the classic introductory text by ( Creutz 1984 1 



We start by writing the partition function of a Eu- 
clidean quantum field theory using the path integral for- 
malism (Dirac 1933 Feynman 1948a|b 1949 Feynman 



and Hibbs 1965 1 as 



and the Euclidean Dirac operator in the case of one 
fermion flavor 



V<^e 



-S(*) 



(1) 



with the action S{^) and <J> genericaUy denoting all fields 
of the theory. For bosonic fields one typically introduces 
periodic boundary conditions, while for fermion fields it 
is natural to introduce antiperiodic boundary conditions 
in time direction (see e.g. appendix A of (Polchinski 



1998)). While this subtlety usually can be ignored, it 



does play some role when choosing the parity of interpo- 
lating operators as discussed in sect. |III.B For a gauge 
theory with fermions one specifically has 



M = jf^Df, + m 



where the covariant derivative is given by 



Df, = df, + igAf, 



(4) 



(5) 



VA., 



X>V2?^e-('^*^'''+^'=) 



(2) 



with the Euclidean gauge action 



1 



Note that in the massless case M is antihermitian and - 
depending on the gauge field configuration - may have ex- 
act zero modes. These zero modes of the operator are re- 
lated to the topology of the underlying gauge field by the 



Sg = ^Ff^^Ff,^ Ff,^ = d^,A^-d^Af,+ig[Af„A^] (3) 



Atiyah-Singer index theorem (Atiyah and Singer 1968) 



Using the rules of Grassmannian integration, one can genericaUy rewrite ^ as 

Z= fvA^ /'l?V^i/;e-('^*^'''+'^'=^ = /"l?A^det(M)e-^<= 



(6) 



Averages that correspond to expectation values of time ordered operators in the operator formalism are genericaUy 
obtained by 



(0> 



1 

z 



VAa / V^jVijOe 



-(4,M4,+Sg) 



(7) 



where O denotes a generic observable composed of gauge and fermion fields. The integration over the fermion fields 
can again be explicitly performed resulting in the replacement of fermion bilinears by propagators on a given gauge 
field background. For a single fermion bilinear this explicitly reads 



V^Wi^{i,Z{y)i^^^{x)) 



-ipMip _ 



det(Af) (M 



CfCi 



{V,x) 



(8) 



while for a general expression it results in the usual combination of all Wick contractions. 



As one can see from (2|8), virtual fermion effects ("sea 
fermions") are contained in the det(M) factor after the 
fermion field integration. Ignoring this det(M) factor 
results in an uncontrolled approximation to QCD, the 



I 

where the a'^^^ are the lattice spacings in direction /i. 
Although more general topologies are possible in princi- 



ple (see e.g. (Jersak et al. 1996)), one usually imposes 



quenched approximation (Marinari et al.\ 1981b). 



Xf^. Here we 



B. QCD regularized on a lattice 

The path integral (IT]) has to be performed over all field 
configurations. In order to make it well-defined, we reg- 
ulate it on a finite spacetime lattice 



^(m) 



with 



e{0,...,iV^-i} (9) 



toroidal boundary conditions x ^ -\- N ^a'^'^^ 
will also specialize to the common isotropic case in which 
the lattice spacing in all directions are equal a^^^' — a. 
The anisotropic case will be separately discussed in sec- 
tion HB 

The fermion field V' is now a Grassmann vector, de- 
fined at the discrete lattice points x — na. We can write 
the naive discretization of the free fermionic continuum 
action as 



with 



and 



^.-li^.-K) 



\^tJ-)xy ~ °x+[i,y 



(11) 



(12) 



and has a naive continuum hmit 

a->-0 - . T „ 1 



C/u. 



1 + iga^F^, - -g^a^Fl, + 0(0^) (19) 



The simplest discretization of the continuum gauge ac- 



tion therefore reads (Wilson 1974) 



Note that this action was obtained by replacing the con- 
tinuum derivative operator d^ with the simple lattice fi- 
nite difference operator T)^. This choice is not unique 
and this non-uniqueness can be exploited to construct 



Sw=PY. (^-\^<Ul{x) + U,,{x)) 



^2 



(20) 



efficient fermion regularizations. Another feature of ( 10 ) 



with (5 — 6/g which has the continuum limit 

a-J-O 1 



5i 



w 



is that it does not describe a single fermion flavor even in 
the continuum limit. The latter is known as the fermion 



d^xTi {F^,{x)F^,{x)) + Oia') (21) 



doubling problem (Karsten and Smit 1981) and will be 
discussed in detail in sect. III.DI In sect. III.E."2] we will 
discuss how the ambiguity in the fermion discretization 
can be utilized to construct numerically efficient lattice 
fermion regularizations. 



We see that ( 20 ) which is known as the Wilson gauge ac- 
tion or plaquette action has discretization errors of O(a^). 
Again, this discretization is not unique and one can uti- 
lize this ambiguity to find gauge actions with higher order 
discretization effects. This will be discussed in detail in 
By combining ^ with ( |11|15| ) and (J20]) 



sect. 



lI.E.l 



The action ([10|) is invariant under a global symmetry and introducing the dimensionless quantities * 



transformation 



* 



,3/2, 



'^, 



il){x) -^ A'ip{x) tp{x) -^ -0(a;)A'f 



(13) 



with A e SU{3) for the case of QCD. This symmetry can 
be promoted to a local one 

?/'(a:) -^ A{x)^^:{x) %i){x) -^ ■4){x)K\x) (14) 

by including a parallel transport Ui^i_{x) to the one-hop 
term V"^ and thus replacing ( 12 ) with 



a'^i'"4>, m = am and D^ 
naive lattice QCD action as 

Nf 

Sn = Sw +22 ^F[mq) 
9=1 



al?,,, we can write the 



Y, {l-\Tr{Ul,{x) + U^,{x)) 



(22) 



{V,),, - U^{x)5, 



+M,1/ 



(15) 



The parallel transport U^{x) is the discretized version 
of the path ordered product of continuum gauge fields 

U^{x) = 5'e'9/x^+'''^^;.^M(^') (16) 

with g being the coupling constant and transforms as 

[/^(a;)^A(a;)f/^(a;)At(x + A) (17) 

under gauge transformations. Note that one could in 
principle choose different paths than the direct one in 



6 

^^§(x)(7^i^^ + m,)*(a:) 

q—l X 



where we have in addition taken the explicit sum over 
Nf fermion flavors q. 

Before we go into the details of the fermion and gauge 
field discretization, let us mention briefly how in prin- 
ciple the cutoff is removed in lattice QCD. As one can 
see from (22), the lattice action exclusively consists of 



(16). As long as the end points remain fixed, the ac- 



tion will be invariant under local transformations ( 14 ) 



This non-uniqueness will play a role when constructing 
efficient fermion discretizations in sect. III. K2] 

In order to construct a kinetic term for the gauge field 
we first note that the trace over a closed loop of paral- 
lel transports is gauge invariant. The simplest of these 
loops, the plaquette, is defined as 

Uf^Ax) = U^,{x)UAx + Mix + iy)Ulix) (18) 



dimensionless quantities. The parameters of the action 
are the fermion masses rriq and the coupling /3. In order 
to remove the cutoff, i.e. to take the limit a — > 0, one 
therefore has to tune these parameters such that on the 
one hand the lattice spacing a goes to zero while on the 
other hand a certain set of dimensionful physical observ- 
ables that are used to define the physical content of the 
theory remain constant. These trajectories in parame- 
ter space of /3 and the m^ along which a set of physical 
observables remains constant as the limit a — >■ is taken 
are called lines of constant physics. A detailed discussion 



of how these can be defined is given in sect. IV 



Along these lines of constant physics it is clear that cor- 
relation lengths in physical units will go to a finite limit 
and therefore will diverge in units of the lattice spacing a. 
In order to possess a continuum limit it is therefore nec- 
essary for a lattice field theory to exhibit a second order 
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FIG. 1 Nonperturbative running coupling constant of 4 flavor 
QCD in units of the QCD scale A from a lattice calculation of 



(Sommer et al. 20101 compared to perturbative calculation. 



Plot reproduced with friendly permission of Rainer Sommer. 



phase transition. Problems can arise if the bare couphng 
constant diverges at a finite cutoff, the so called Landau 
pole ( Landau 1955 ) . In that case the only line of con- 



stant physics that does not show a divergence at finite 
cutoff is the one with vanishing coupling, i.e. the trivial 
theory. In order for theories with a Landau pole problem 
to have a non vanishing renormalized coupling (i.e. to be 
nontrivial) one must retain a finite cutoff which prevents 
one from taking the continuum limit. Such theories can, 
however, still serve as effective theories. Consequences 
for the lattice formulation of this class of theories are 



discussed e.g. by ( Arnold et al. 2003 Espriu and Taglia 



cozzo[ [2003} [Gockeler et al. 1998a Kogut and Strouthos 



2Q05D 

Due to asymptotic freedom (Gross and Wilczek 1973 



Politzer 1973) however, no such problems are expected 
to arise in lattice QCD. The perturbative expectation 
of the vanishing of the QCD coupling constant at large 
scales has been confirmed by nonperturbative lattice cal- 



culations in various settings (Bode et al. 2001 Bowler 



et al. 1986 Delia Mor te et al] [2003 
Luscher et al.\ ,1994^ |Sommer et al 



[Gupta et al. 



1988 



2010 Tekin et al. 



2010[ ). The most recent result for Nf 
ted in fig. [l] 



4 QCD is plot- 



C. Numerical evaluation of the Path Integral 

Before we return to the task of constructing a lattice 
regularization of QCD, we need to discuss some basics 
of the numerical evaluation of the path integral ((7|. In 
terms of dimensionless lattice quantities, ([7| can be writ- 



ten as 



(0) 



1 

z 



VUa / X>*X>*Oe- 



-(*M(C/)*-|-Sg((7)) 



(23) 



with a generic gauge action Sg{U) and fermion operator 
M{U), which in general depends on the gauge field U. 
We would now like to perform the integration over both 
the fermion fields ^, \1/ and the gauge filed U. Using 
the relation ([8]) and its generalization for fermion mul- 
tilinears, we can explicitly perform the integration over 
fermion fields. With the understanding that we have to 
replace fermion multilinears by the sum over all Wick 



contractions, we may thus cast (23 1 into the form 
1 
Z 



(O) 



Yl dU^{x) det M([/)Cle-^« 



(U) 



(24) 



Generally and for QCD in particular, it is not possible to 
perform the remaining integration over the gauge fields 
U^ in closed form. While both strong (Wilson 1974) 
and weak (see e.g.(Capitani 2003)) couphng expansions 
are possible, neither allows for a detailed quantitative 
understanding of the nonperturbative dynamics of the 
system. 

From a numerical perspective, ([24[) is a high dimen- 



sional integral over 4 x J| N^ copies of the gauge group, 
which in the case of QCD is S'C/(3). The only category 
of numerical methods that are suitable to perform such a 
high dimensional integral are stochastic or Monte Carlo 
(MC) methods in which the space to be integrated over 
is randomly sampled, i.e., where observables are aver- 
aged over on randomly drawn gauge field configurations 
Ufi{x). As Monte Carlo integration is stochastic in na- 
ture, there is always a statistical error associated with it. 
This error has to be estimated in a lattice calculation, 
which is typically achieved via the jackknife or bootstrap 



methods (see e.g. (Press et al. 

A straight Monte Carlo integration of ([24[) however 



2007)) 



in which the gauge configurations C/^ are randomly pro- 
duced with equal weight is extremely inefficient. For in- 
teresting parameter choices, all but a very small sub- 
set of relatively smooth configurations are exponentially 
suppressed by the exponent of the gauge action Sg{U) 
and the fermion determinant det M{U). In order to 
circumvent this problem, one can produce gauge field 
configurations with a probability that is proportional to 
det M([/) X e"'^'^'-'^-' and compute the expectation value 
of an observable as an unweighted average over these con- 
figurations. This technique, known as importance sam- 
pling, requires an algorithm that produces gauge field 
configurations with the proper weight. Typically this is 
achieved via a Markov chain process, where a time series 
of gauge field configurations is produced in which the n"^ 
configuration Up^ depends on the previous one C//t" 

A Markov chain is characterized by the transition prob- 
ability 



p(m,n)=P(C/™|[/") 



(25) 



where the J7* are all possible gauge configurations and 
the conditional probability P{U"^\U") is understood in 
the sense that it denotes the probability of the system 
to go over from a configuration U" to U™ in one time 
step. The transition probability p acts on the space of all 
gauge configurations. It fulfills the two basic relations 



Vm, n : p{m, n) > / dmp{m, n) = 1 



(26) 



where J dm denotes the integration over all possible 
gauge field configurations. If in addition the transition 



probability ( 25 ) fulfills the detailed balance condition 
Vm, n : p{n, m)pm = p(m, n)p„ 



(27) 



with the desired equilibrium distribution p„ — 
det M(C/") X e~^'^'^'^"yz, then one can easily show the 
following two properties: 

1. The transition probability maps the equilibrium 
distribution onto itself 



Prn = / dnp{m,n)p„ 



(28) 



2. Defining a distance d{w^v) = J dn\wn — w„| in the 
space of probability distributions, the application 
of the transition probability moves every probabil- 
ity distribution closer to the equilibrium distribu- 
tion 



d{pw,p) < d{w,p) 



(29) 



If in addition to (28) and (29) the system is ergodic, 
i.e. if any configuration C/' may be reached from any 
other configuration U^ with non- vanishing probability in 
a finite number of time steps, then it is guaranteed that 
starting from an arbitrary initial probability distribution 
we end up with the desired equilibrium distribution p. 

The time until the equilibrium distribution p is reached 
(in the sense that no statistically relevant drift towards 
the equilibrium expectation value can be seen in any 
monitored observable) is usually called thermalization 
phase and its shortness is an important quality crite- 
rion of an algorithm. Once the system is thermalized, 
i.e. the equilibrium distribution has been reached, it is 
advantageous if consecutive configurations have as little 
correlation as possible. In order to have a quantitative 
handle, it is customary to monitor the autocorrelation 
time of certain observables within a Markov chain. 

In the case of a pure gauge theory or the quenched 
approximation where the fermion determinant factor 
det M(U) is missing and the weight factor is proportional 
to the exponent of the gauge action e~^'^^^\ the update 
algorithms that produce the next element in the Markov 
chain usually exploit the locality of the gauge action 
Sg{U). As pioneered by ( Creutz et al. 1979a|b ), one can 



pick a certain gauge link Uf^{x) from the current gauge 
configuration U and produce a suggested new gauge con- 
figuration U' by multiplying [/^ {x) with an element of the 
gauge group. Since the gauge action Sg{U) is a sum of lo- 
cal terms, the change in the action SS = Sg{U') — Sg{U) 
is readily evaluated by recomputing those few terms that 
contain the flipped gauge link. One can then perform 
a Metropolis (Metropohs et al. 1953) step, i.e. accept 



the gauge configuration U' as the next gauge configu- 
ration in the Markov chain with probability e~ if the 
action has increased SS > or with probability 1 other- 
wise. It is readily seen that this algorithm satisfies the 



detailed balance condition (27). Another frequently used 



local update algorithms for pure gauge theories is the 
heatbath ( |Creutz[[T980b} [Kennedy and Pendletollll985|) 
Supplemented by overrelaxation steps ( Adler 1981 1988 



Brown and Woch[|1987 Creutz 1987 Fodor and Jansen 



1994), these algorithms are still the state of the art for 



pure gauge theories. 

Due to the nonlocal nature of the fermion determinant 
det M{U) an update in a theory with dynamical fermions 
is substantially more complex and computationally de- 
manding. For a lattice with N = J| N^ sites, M{U) is 
typically a (12 x N)"^ matrijFland therefore a direct com- 
putation of detAf([/) is prohibitively expensive for even 
moderately sized lattices. Although alternative sugges- 
tions have been made (Berg and Forster 
T982; Luscher 



et al. 19811 Kuti 



1994 



1981 



Montvay 



Polonyi and Wyldl ,1983) |Scalapino and Sugar 



Fucito 



1984 



1981 



Slavnov 1996), one usually proceeds by introducing a 



bosonic (complex scalar) pseudofermion field 4> ( Wein- 
garten and Petcher 1981[ ) . The fermion determinant may 
thus be written as 



det M{U) 



P$^P$e-**^^(^)"'* 



(30) 



The catch here is of course the appearance of the inverse 



fermion matrix M{U)~^ in (30) which is again a nonlocal 
object. In addition, the kernel operator M{U)~'^ has to 
exist (i.e. the matrix M(U) needs to be invertible) and 
be positive definite hermitian in order to ensure the con- 
vergence of all gaussian integrals over the pseudofermion 



field in (30). From (22) we see however, that the naive 



fermion operator is not hermitian and neither will be the 
fermion operators we will construct later on. As long as 
det M is real and positive definite however, one may use 
the identity det M = ^/det{M^M) to rewrite an arbi- 
trary power of the fermion determinant as 



det Miuf 



P$^2?$e-*H^^^(^)^^(^))" 



(31) 



^ Note that in the staggered fermion formulation the size of the 
matrix is reduced to (3 X A'^)^. For a detailed discussion see 
sect.jlLDl] 



The path integral can now be formulated in terms of 
bosonic variables only with an additional term in the 
action 



Sf = $^ {m\U)M{U)) 



-a/2 



$ 



(32) 



Note that for actions where M\U)M{U) does not couple 
even and odd lattice sites one may choose to keep the 
pseudofermion fields <i>e on even lattice sites only and 
thereby obtain 



det A/([/)" = / P$t2?$,e-*'(*^'(^)*^(^))" 



In order to efhciently integrate the system of pseud- 



ofermions and gauge fields we follow (Batrouni et al. 



1985 Callaway and Rahman 1982, 1983; Duane 1985 



Duane et aL||1987[|Duane and Kogutnl985Hl986iiPolonyi| 



and Wyld 1983 1 and reinterpret the total action of the 



system 



Sg + ^^ {m\U)M{U)) 



-a/2 



$ 



as the potential part of a fictitious Hamiltonian 



^-V 



5(0) 



(34) 



(35) 



with conjugate momenta tt where collectively denotes 
all pseudofermion and gauge fields. One can then pro- 
ceed to choose some initial momenta and integrate the 
canonical equations of motion 



dS 



(36) 



numerically in a fictitious time r along a "classical" path. 
The classical partition function corresponding to the set 
of all such classical trajectories is given by 



X>7rX></«e 



-H 



Vne 



V<j>e- 



(37) 



As the gaussian integration over the momenta tt only 



gives an irrelevant prefactor, ( 37 1 does reproduce the cor- 



rect probability distribution in the original theory. As- 
suming ergodicity, one can obtain the correct distribu- 



tion of classical paths ( 37 1 by periodically refreshing the 



momenta tt with a random value from a gaussian distri- 
bution. The expectation value of an observable can thus 
be obtained by averaging it along all classical trajectories 
in the update chain. The inexact nature of the numerical 
integration introduces a systematic error, which however 
can be corrected by a final Monte Carlo accept/reject 
step of the complete trajectory (see (Duane et al. 1987)). 
This algorithm is known as hybrid Monte Carlo (HMC). 



The Hamiltonian in ( 35 1 is readily constructed in the 
case where a in (31) or (33 1 is a positive integer. For a 



general fractional power a one can resort to a polynomial 



(33) (de Forcrand and Takaishi 1997 Frezzotti and Jansen 



1997) or rational (Clark and Kennedy, 2004 Clark et al. 



2005 ) approximation of the desired fractional power of 
W{U)M{U). These versions of the HMC algorithm are 
known as polynomial HMC (PHMC) and rational HMC 
(RHMC) respectively. 



In order to integrate the equations of motion ( 36 ) nu- 



merically, one has to compute the derivative of the action 



( 34 ) with respect to the gauge field 



dS _ dSa 



$ 



d{M\U)M{U)Y 



$ (38) 



This derivative is commonly known as the force term. 
The computationally expensive part of ( |38[ ) is the second 
term, the fermionic force term. In the special case a = 1 
and using the shorthand notation M{U) ^ M^{U)M{U) 
it may be written as 






-^^M{U)-'^^^^MiUr'<^ (39) 
dUf,(x) 



We see that a single inversion of A4{U) on a pseudo- 



fermion vector is required to compute ( 39 ) . For general 



fractional exponents a one may introduce a rational ap- 
proximation 



with 



r(7W(t/))~X(C/)-" 



rix) = Yl 



x + Pi 
In this case, the fermion force may be written as 



(40) 



(41) 



$V(X(C/))$ = -^$^ 



dUf,{x) 



)-$ 



(42) 



which can be computed using one single multilinear ma- tor. 



trix inversion (de Forcrand 1996 Frommer et al. 1995 



Glassner et al. 



1996 Jegerlehner 1996 ) on a single vec 



An alternative integration scheme, the hybrid molecu- 



lar dynamics (HMD) R-algorithm, was proposed by (Got 



tlieb et al. 1987b). Although it has seen considerable use D. The fermion doubling problem and its solutions 



in the past, it has largely been replaced by the RHMC al- 
gorithm. It is a pure molecular dynamics algorithm that, 
in contrast to pseudofermion algorithms, does not allow 
for a final MC step to correct for finite step size errors ac- 
cumulated along the integration trajectory. Due to this 
feature, detailed balance is fulfilled by the R-algorithm 
only in the limit of a vanishing step size in contrast to 
HMC-type algorithms that are exact even at finite step 
size. 

An efficient HMC algorithm has to simultaneously sat- 
isfy two criteria: On the one hand, the acceptance rate 
should be high (one typically aims for '^ 80 — 90%) and 
on the other hand the autocorrelation time should be 
small. The autocorrelation between successive configura- 
tions can be decreased by a longer integration trajectory 
separating them. This however leads to larger numeri- 
cal integration errors and consequently to a lower accep- 
tance rate. A trivial remedy consists of decreasing the 
time step in the numerical integration, which however is 
computationally expensive because the fermion force has 
to be computed more often. It is therefore advantageous 



to use higher order integration schemes ( Omelyan et al. 
2002a|b[ [2003) [Takaishi and de Forcrand[ |2006[ ) that al 



low a larger time step in the numerical integration while 
keeping the acceptance rate high. 

Another method for speeding up HMC type algorithms 
consists of introducing different time steps for pseud- 
ofermions and gauge fields (Sexton and Weingarten 



1992). Splitting off the UV modes of the spectrum 



and Jansen 



2003 



by mass preconditioning ( Hasenbusch 2001 Hasenbusch 



or via domain decomposition ( Luscher 



2003[ |2004[ |2005[ ) and integrating IR and UV part with 
different time steps leads to a substantial additional 
speedup. This speedup is especially large if combined 
with the suppression of UV modes and other improve- 
ments of the fermion regularization that will be discussed 
in sect. lII.E.2l 

As noted in sect. |II.A[ ignoring the effects of the 
fermion determinant results in the quenched approxima- 
tion. Since it bypasses the most computationally de- 
manding part of the ensemble generation, it was exten- 
sively used in the early years of lattice QCD and is still 
useful for certain conceptual studies. Although it is an 
uncontrolled approximation, it may be justified by not- 
ing that it becomes exact in the large N,, limit. Fur- 
thermore, by choosing the proper scale setting observ- 



able (see sect. IV) a large part of the dynamical fermion 
corrections might cancel and effectively be absorbed into 
a redefinition of the coupling constant |j 



We now return to the free, naive fermion action (22) 



* (7m ^P + ™) * 
The fermion operator reads 

M = 7^13^ -I- m 

which in Fourier space becomes 



M{p) = -^1^, sin(ap^) 



(43) 



(44) 



(45) 



The momentum space propagator is consequently given 
by 



D{p)=M-'{p) 



iEu7Msin(aPp 



(iEpSin(app 



(46) 



has 



which in addition to the physical pole at p^ — - 
15 additional poles located at the edges of the Brillouin 
zone. The poles are located at {p — 11)^ = —rn^, where 
n is any of the 16 four-momenta 



n= (po,Pi,P2,P3) with p^ e {0, 7r/a} 



(47) 



This rather fundamental obstacle of putting fermion 
fields on the lattice is known as the doubling problem. 
Physically, we can trace this problem back to the well 
known axial anomaly of a continuum theory. In the mass- 
less limit, a classical fermionic theory is invariant under 
the chiral transformation 



*(a;) -> *'(a;) = e''^'^^*(x) 



(48) 



As demonstrated by ( Adler 1969 Bell and Jackiw 1969 1 , 



the conservation of the corresponding Noether current, 
the axial vector current, is destroyed by quantum fluc- 
tuations. In a lattice regulated theory however the ex- 
istence of a classical symmetry implies a conserved cur- 
rent. The anomaly of the physical fermion axial vector 
current is canceled by the anomaly of unphysical doublers 



1981 



as demonstrated by (Karsten and Smit 

It was later shown by (Nielsen and Ninomiya 



1981a|b|c), that no lattice fermion regularization exists 



that fulfills all of the following conditions at the same 
time 

i Absence of doubler fermions 



* Another attempt to ju stify the u se of the quenched approxima- 
tion was made by ( fAnthony et al.\ \l982 . Duffy et al. . 1983 ). They 
gested to extrapolate to a positive number of quark flavors 



sugge 



by computing observables in the quenched approximation and at 
an effecti ve negative number of quark flavors. 
5 See also JChodos and Healyl |l977l iKerlerl |1981 1 



(a) Ferinionic case 

. . < nonlocality 

, doubler fermion 




(b) Bosonic case 



FIG. 2 In the fermionic case |(a)| periodicity of the lat- 
tice momentum P^ requires a second root (doubler) or jump 
within the Brillouin zone (nonlocality). In the bosonic case 
(b) only the squared lattice momentum is required to be pe- 
riodic which can be fulfilled without a jump or an additional 
root. 



ii Continuum chiral symmetry in the massless case 

ill Locality in a sense that M{x,y) — > vanishes expo- 
nentially as X — y ^>- oo 

iv Correct continuum limit 

This result can be understood by noting that a general 
lattice fermion operator M which anticommutes with 75 
in the m = case can be written as 

M{p)=m + iJ2 7m ^M («P) + J2 7a'75-Rm (ap) (49) 

The requirement that it reproduces the correct contin- 
uum theory implies, that for small a the lattice momen- 
tum P^ goes over into the continuum momentum p^ while 



which only depends on the discretized momenta squares 



(^\-:i(i 



iaPfj.) 



(52) 



These are naturally periodic with a period of 27r/a as 
displayed in fig. [2b] 

One therefore has to give up on any one of the above 
requirements for lattice regularizations of fcrmions. Ob- 
viously, one can not give up the requirement ( pv| ) of a 
correct continuum limit. Giving up the locality require- 
ment mw on the other hand has been suggested, among 




others, by ( Drell et al. 1976 ) who proposed 



Pf^^Pf^ 



by (Rebbi 1987), who suggested 



2V sin^ 

P^ = sm ap^ -^= — -. 

a 2^^ sm ap„ 



(53) 



(54) 



and by (Gross et al. 1987) whose construction involves 



non-symmetric difference operators and contains non- 
renormalizable terms in the continuum limit. However, 
all of these approaches turned out to be problematic and 
have been abandoned. 

Among the remaining two options, we will first dis- 
cuss latice fermion regularizations that give up on the 
requirement m and therefore describe more than a sin- 



gle flavor in the continuum limit. Among these, (Borici 



2008 Creutz 



2008 Karsten 1981 Wilczek 1987) have 



suggested different implementations of minimally dou- 
bled fermions, i.e. lattice fermions which have one single 
doubler only. As this single doubler has to be placed 
somewhere within the Brillouin zone, all of these for- 
mulations share the characteristic that in Fourier space 
there is a distinguished direction, namely the direction 
from the physical particles pole to the pole of the single 
doubler fermion. Therefore, a number of discrete lattice 



symmetries are broken ( Bedaque et al. 2008 ) resulting in 



a more complicated renormalization pattern and, gener- 
ically, in a fine-tuning of the parameters of the action 



D V n Ajj'-t- 11 n • • J- ■ J- i- (Capitani et al. 2010). Currently fundamental proper- 
Rfj, — >■ 0. Additionally, -Pp is periodic m every direction 1 — i- _- r I 1 . -^ . . . 



with period 27r/a. As shown in fig. [2] these restrictions 
on P^ imply that either it has a second root in the first 
Brillouin zone, which gives an additional pole in the prop- 
agator, i.e. a doubler fermion, or that it has at least one 
discontinuity, which makes the fermion operator M{x,y) 
nonlocal. Note that in comparison the discretized version 
of the continuum action of a scalar field 



Sr. 



= 2^^ {d,d^ - ^^) </> 



(50) 



in momentum space reads 



ties of minimally doubled fermions are still being clarified 
and applications to hadron spectroscopy or other phe- 
nomenologically relevant computations are not available 
in the literature yet. 



1. Staggered fermions 

A less minimal but more symmetric way of putting 
doublers on the lattice is given by the staggered fermion 
formulation that was developed in a series of papers by 



Ss — — 



1 






p-2 



(Banks et al. 1976 Kogut and Susskind 1975 Susskind 



(51) 1977). Staggered fermions are obtained from the naive 



fermion action ( 43 1 by noting a fourfold exact degeneracy 
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(in the interacting theory) that can be exposed by a spin- 
diagonahzation 



^{x)^r{x)x{x) ^{x)^x{^)THx) 



with 



r(x) = n 



7/ 



(55) 



(56) 



In terms of x ^-nd Xi (43 1 can be written as 



x{Vt.D^ + m)x Vt.{x)^{-1)^''>^"- (57) 

where 'ri^{x) is a pure phase factor making exphcit the 
decouphng of the 4 spin components of x- Defining x on 
a single component only, we have reduced the fermion 
content of the theory by a factor of 4, from 16 four- 
component spinors to 16 single component modes that 
are still symmetrically distributed over the Brillouin zone 



at the momenta 11 given in ( 47 1 



The staggered fermion operator 
M = Tj^,D^ + r 



(58) 



is antihermitian in the massless case, i.e. its eigenvalues 
are restricted to the imaginary axis for m — 0. Therefore 
any finite mass m > provides an IR cutoff and the 
operator is invertible. In addition, the massless staggered 
fermion operator does preserve a remnant of the chiral 



symmetry of the naive fermion operator ( 48 1 



X{x) -^ x(x) = e^^^xix) 
X{x) -> x{x) =x{x)e'^'' 



(59) 



where e = (—1)^'' ^^ . For the staggered fermion operator 
( 58 1 this implies that 



eM = M^e 



(60) 



is hermitian. Therefore, the eigenvalues of M come in 
complex conjugate pairs and det M is real and positive. 
Furthermore, we see that 



MHd = r]^Dlr],D, + \m\'' 



(61) 



does not couple odd and even sites so that the pseudo- 



fermion representation ( 33 1 may be used for the fermion 
determinant. 



As demonstrated by (Daniel and Kieu 



Doel and Smit 1983[ Gliozzi 1982 [Golterman and Smit 
1984b Kluberg-Stern et al. 1983| [Sharatchandra et al. 



1986 



van den 



1981), these 16 spinor components may again be inter- 



preted as 4 fermion flavors (also referred to as tastes 
in the literature), each one described by a 4-component 
spinor. Following ( Golterman[ 1986), we split the lattice 
coordinate x = y + h into a piece y that describes the 



origin of the the elementary 2^ hypercube that x is lo- 
cated in and an offset h which describes the location of 
X within this hypercube. We construct 

X{y) = \Y.^{h)U{y + h,y)x{y + h) (62) 



where ft.^ e {0, a} and U{y+h, y) is any parallel transport 
from y + h to y. Interpreting X{y) as a 16 component 
vector, we can define an arbitrary fermion bilinear with 
spin structure 7=, and flavor structure 7/ as 



(63) 



X{^,^^f)X = TT[X-f,X-f} 



and the staggered fermion action (57 1 can be written as 



X 



(<■ 



1)™+(7m®1)^m + (75®Ca.^5)C7^U (64) 



with the first and second derivative operators on the 
coarse lattice 



D„ 



where 



\i^'-n 



Ca 



1 



K 



21 + V;n (65) 



{-^ 



Uf_,{x + fl)Uf,{x)Sx+2fi,v 



(66) 



Note that the third term in ( 64 ) that implies a mixing 



of the four remaining flavors (tastes) is an artefact of 
the taste assignment ( 62 1 which does not respect the full 
set of symmetries of the staggered action. In fact, for 



the noninteracting case ( Adams 2005 ) has found a taste 



assignment that is diagonal in the tastes and local. 

A massless four flavor continuum theory does have a 
classical f/(4) chiral symmetry which gets reduced to 
5L/(4) by the anomaly. This implies a 15-plet of mass- 
less pseudoscalar Goldstone particles (pions) and an ad- 
ditional massive one (77'). As we have seen in (59), stag- 
gered fermions retain a U{1) subgroup of this symme- 
try. In the spin-flavor basis, the remnant staggered chiral 
symmetry reads 



Xiy) 



— pitPill, 



X'{y) 

X'{y)=X(y)e 



^'^^X{y) 

10(758155) 



(67) 



This symmetry is spontaneously broken implying a single 
Goldstone particle, the pseudoscalar in taste space, that 
is exactly massless. 

In order to obtain one resp. two flavors in the func- 



tional integral (24 1, it is customary to take the quartic 



resp. square root of the 4- flavor staggered functional de- 
terminant det M. This procedure, commonly referred to 



as rooting, was introduced in (Marinari et al. 1981b I in 



the context of the Schwinger model. On a technical level 
this is realized in a pseudofermion based algorithm by a 
fractional power a = 1/4 resp. a = 1/2 in (33 1. 
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On a more fundamental level, the validity of rooted 
staggered fermions relies upon the assumption that there 
exists a local lattice fermion operator that squared or 
to the fourth power has the same functional determi- 
nant detM as the 4-flavor staggered operator. In fact, 



( Adams 2005 1 demonstrated that in the free case such an 
operator can be found. In the interacting case however. 



the situation is more complex. As (Durr and Hoelbling 



2005b ) have demonstrated, rooted staggered fermions are 



in the wrong universality class in the strictly massless 
case TO = implying that the chiral m — )■ limit does 
not commute with removing the cutoff a — > 0, a result al- 
ready anticipated by (Smit and Vink 1987). This obser- 



vation reflects the fact that the staggered chiral symme- 
try ( 59|67 ) is not anomalously broken. In fact, the stag- 
gered chiral symmetry is retained even in the single flavor 
rooted theory implying an exactly massless 77' for to = 



and as demonstrated by (Bernard 2006 Bernard et al. 
2007a|b Prelovsek 2006 ) nonunitarity of the rooted the- 
ory at finite cutoff. On the other hand, there is some 
numerical indication that out of the chiral limit rooted 
staggered fermions do indeed have correct continuum be- 



havior for many observables ( Aubin et al 



eTall |2010a1bl [Davies eFar] |2005[ |2004[ |Durr and Hoel 



bling |2004 2005b Durr et al. 2004 FoUana et al 2008). 



More formally, (Bernard et al. 20061 have shown by a 



2004 Bazavov 



symmetry argument that rooted staggered fermions can 
not be described by a local operator. However, analyti- 



cal calculations ( Bernard ei aL 2006 2008a Giedt,2007 



Shamir 2007 ) indicate that the nonlocal terms vanish in 



the continuum limit and consequently that rooted stag- 
gered fermions do have the correct continuum limit as 
long as the proper order of limits is observed. Impli- 
cations of the delicate nature of the staggered fermion 
continuum limit were also extensively discussed in the 



2007a 



mann 



literature ( Bernard 20051 Bernard et al 2007b| Creutz 



b) Durr and Hoelbling" '2 006| Hasenfratz and Hoff-| 
2006). For recent reviews see e.g. (Durr, 2006 



Sharpe 2006) 



The representation of the spin-taste structure by differ- 
ent points within an elementary hypercube ( 56|62 ) and 
the taste breaking term in the action ( 64 1 imply some 



additional complications for extracting hadron masses 
with staggered fermions that will further be discussed 
in sect. IIIII 



2. Wilson fermions 

We now turn our attention to lattice fermion formu- 
lations that fully lift the naive flavor degeneracy and 
are able to naturally describe a single flavor theory in 
the continuum limit. It was first realized by (Wilson 



1975) that the fermion doubling problem can be solved 



by adding a laplacian term to the naive fermion operator 



(43) 



Sw^^h^.D^ 



n]^ 



where 



D 



E^^ 



C^^V^-2 






(68) 



(69) 



with the parallel transport V^ defined in (15) and the 



Wilson parameter r that is usually set to 1. The addi- 
tional term in the action, the so-called Wilson term, may 
be interpreted as a momentum dependent mass term. 
Note that in contrast to naive and staggered fermions, 
the Wilson fermion operator is generally not normal due 
to the additional laplacian term, although it still is in the 
free case. The free Wilson operator 



Mw = 7m^p 
in momentum space reads 



-D 



(70) 



Mw{p) = - E ^^ sin(ap^) + to ^ (cos(ap^) - 1) 



(71) 



r 



Comparing ( 71 ) to the naive operator (45 1 we see that the 
additional term ^ ^ (cos(ap^) — 1) does vanish as 0{a) 
for any fixed physical momentum p. On the other hand, 
for a fixed lattice mom,entum, ap the additional term gives 
a contribution that is divergent as 0{l/a) except for p — 
0. In particular, all doubler modes with n momentum 
components 7r/a do receive an additional mass of Irnja 
thus effectively removing them from the spectrum in the 
continuum limit. 



Since the additional laplacian term in ( 68 ) does not 
anticommute with 75, the exact chiral symmetry of 



naive fermions ( 48 ) is broken as required by the Nielsen- 
Ninomiya theorem ( Nielsen and Ninomiya 1981a|b|c ) . In 
fact, the laplacian term commutes with 75 due to its triv- 
ial spin structure. Consequently, the Wilson operator 
obeys the relation 



M, 



w 



I5MWI5 



(72) 
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which is known as 75-herniiticity. It impHes that the 
operator 



75MH/ = {-i^Mw) 



(73) 



is hermitian and the eigenvalues of Mw are either real 
or come in complex conjugate pairs. Consequently, 
det(Mtv) is real. In order for the pseudofermion rep- 
resentation (31| to be well defined, de\.{Mw) needs to 
be positive definite in addition, which is guaranteed if 
TO > 0. However, due to the breaking of chiral symmetry 
the fermion mass is not protected against additive renor- 
malization. The bare fermion mass receives corrections 
that are divergent in the continuum limit and needs to be 
renormalized. Due to this additive renormalization, the 
bare fermion mass corresponding to a physically inter- 
esting renormalized mass often turns out to be negative. 
While pairs of complex conjugate eigenvalues still give a 
positive contribution to det(Mvi/) in this case, the real 
eigenmodes only do so if the number of negative ones 
is even. Furthermore, even eigenmodes that are positive 
but very small pose serious problems for matrix inverters. 
In the case of an exact zero eigenmode, it is not possible 
to define the fermion determinant in terms of the pseudo- 



fermion fields according to (31). Hitting such a configu- 



ration within numerical precision will result in a failure 



of the matrix inversion to properly converge (Bardeen 



et al., 1998). These configurations are known as "excep- 



tional" and the appearance of even a single exceptional 
configuration in a Markov chain indicates that one is not 
able to properly sample a relevant region in configuration 
space. Ensembles exhibiting an exceptional configuration 
therefore have to be discarded. 

In practice, exceptional configurations therefore set a 
lower limit to the masses one can reach with Wilson-type 
fermions. One has to make sure that all eigenmodes of 
the fermion matrix are sufficiently separated from zero. 
While these restrictions were initially very severe, they 
do not present a substantial obstacle for current state 
of the art lattice calculations. The use of large physical 
volumes, small lattice spacings, improved gauge actions 



(see sect. lI.E.l) and smeared link fermion actions (see 
sect. II. E. 2 ) all reduce the probability of exceptional con- 



figurations appearing in a simulation. 

Due to the strong correlation between the condition 
number of the fermion matrix and the iteration count of 
the inverter and because the relative fluctuations of the 
largest eigenvalue are small, one can use the distribution 
of the inverse iteration count instead of the distribution 
of the lowest eigenmode. A tail of this distribution that 
extends towards the origin is a clear and direct indica- 
tion of problems with exceptional configurations while 
a clear separation from demonstrates the absence of 
exceptional configurations and positivity of the fermion 
determinant. Such a distribution is plotted in fig. [S] for a 



recent study with light Wilson- type fermions ( Durr et al. 



2011b) 



Inverse iteration count (1000/N„) 



|3=3.31, M„-135MeV 



(3=3.5, M„-130MeV 



=3.61,M„-120MeV 



X 



|3=3.7, M„-180MeV 
|3=3.8, m!,-220 MeV 




0.04 0.08 0.12 0.04 0.08 0.12 

FIG. 3 Inverse iteration count of the fermion matrix inverter 
for Wilson-type fermions and a number of different ensembles 
with bare couplings /? and approximate pion masses A/^. The 
lower tail of the distributions shows a clear separation from 
indicating the absence of exceptional configurations in the 
ensembles. From (Durr et al. 2011b I 



The dominance of low modes in the computational cost 
of inverting a Wilson-type Dirac operator has lead to 
efforts of preconditioning the inversion by removing or 
effectively projecting out a relatively small number of 
low modes. These techniques are generally known as 



deflation methods (Darnell et al. 



Schaeferl [2004) [de 



200m 



Luscher 



[2007 



2008 DeGrand and 



tarcrand 1996| Giusti et al.\ 



2004 



Neff et al. 2001 Stathopoulos anc 



Orginos 2007 1 . They can lead to a huge decrease in the 



cost of computing propagators on gauge configurations, 
especially in circumstances where one needs to compute 
many propagators on the same gauge configuration. In 
this case, one can perform a rather expensive eigenmode 
projection step since it only has to be performed once for 
each gauge configuration. For a further discussion see 
also sect. IIII.BI 

On a more fundamental level, real modes of the Wilson 
operator are related to the topological charge of the gauge 



configuration by the index theorem (Atiyah and Singer 
1968). In the continuum limit, real modes become ex- 



actly degenerate zero modes that are tied to the gauge 



field topology ( Bardeen et al. 



T9981 [Hasenfratz et aLI|1998||Setoodeh et a/.||1988[|Smit| 



and Vink 1987 Vink 



1988) 



1998 Gattringer and Hip 



Recently it has been proposed to construct a Wilson- 
like operator that instead of lifting the flavor degeneracy 



of the naive operator ( 44 ) does lift the taste degeneracy 



of the staggered operator (^58m Adams 2011 de Forcrand 



et al.\ |2010[ |Hoelbling[|201iP (for earlier work in this di 



rection see also (Becher and Joos 1982 Gockeler 



Golterman 1986 Golterman and Smit[ |1984a|b 



1983 Mitra and Weisz 19831). Conceptual aspects of 



1984 



Mitra 



this formulation are still being studied and no applica- 
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tion to hadron spectroscopy or other phenomenologically 
relevant computations are available in the literature yet. 



3. Twisted mass fermions 



Twisted mass fermions (Frezzotti et al. 2001) are a 



variant of the Wilson fermion formulation that has re- 
cently gained attention. The basic idea is to perform a 
chiral rotation, that is not affected by an anomaly, on 
the mass term. Because the transformation has to be 
anomaly free, the number of flavors to be chirally ro- 
tated has to be even. In the simplest case of two flavors 
the mass term reads 



where 



j^g«Q!75T3 ^ T^ J__ j^^^^^ 



tan(Q;) — — m = m + fi 

m 



(74) 



(75) 



and T3 is the diagonal Pauli matrix in flavor space. Due 
to the opposite twist angles between the two flavors, the 
anomaly cancels and the chiral rotation of the mass term 



( 74 1 may be absorbed into a chiral rotation of the fermion 



fields 



* -^ *e 



ia/275r3 



qi -^ g«"/275T3v[; 



(76) 



provided that the massless part of the fermion operator 
is invariant under the chiral transformation (76 1. For 
Wilson fermions however, chiral symmetry is explicitly 
broken. Replacing the standard mass term in ( 70 ) with 



a twisted mass term of the form ( 74 ) therefore results in 



a different theory where the two flavor fermion matrix is 
given by 



Mtr. 



7p 



D„ 



+ -D + TO + i^iTsTs 



This represents the twisted mass fermion matrix in the 
so-called twisted basis. The basis is called twisted be- 



cause ( 77 1 does describe the physically uninteresting case 
of a complex mass term. In order to obtain physically in- 
teresting predictions for a theory with a real mass term 
from (77 1, the chiral rotation in the mass term has to 



be supplemented by an equivalent transformation of the 
fermion fields (76). In the new basis, (77) describes the 
physically interesting case of real mass fermions. This 
rotated basis of the fermion fields is therefore usually re- 
ferred to as the physical basis. 

Wilson fermions with a twisted mass term do not obey 
standard time reversal and parity transformation sym- 
metries but modified versions thereof. However, stan- 
dard CPT symmetry is fulfilled and the behavior with 
regard to chiral and flavor symmetry is the same as for 
standard Wilson fermions, i.e. chiral symmetry is broken 
while flavor symmetry is exactly conserved in the twisted 



basis. In the physical basis ( 76 ) however a subset of the 



flavor and axial symmetries get transformed into each 
other. Consequently, flavor symmetry is broken while 
part of the chiral symmetry is restored at maximal twist 
a = 7r/2. This implies that at maximal twist on the one 



hand there are isospin breaking cutoff effects ( Bar 2010 
Scorzato[ [2004 ) while on the other hand the cutoff terms 



are generally of 0{a^) ( Aoki and Bar 



2004 



Frezzotti and 



Rossi 2004a). Note, that the bare mass is not protected 



against additive renormalization and therefore the mix- 
ing angle a gets renormalized. In order to achieve maxi- 
mal renormalized twist, the bare mass needs to be tuned. 
This tuning is routinely done as part of any twisted mass 
calculation (see e.g. (Baron et al. 2010b)). 



Introducing a pair of non-degenerate quarks is usu- 
ally done with the help of an additional mass term that 
carries a nondiagonal flavor structure ti resulting in a 
non-degenerate 2-flavor fermion operator of the form 



Mt, 



ltJ.Df, + -D + m + ifJ,'^5T3 + eri (78) 



(Chiarappa et al. 2007 Frezzotti and Rossi 2004b). An 



alternative suggestion has been proposed by (|Pena et al. 



2004) 



From (77 1 it follows, that the two flavor operator has 
a determinant which is bounded from below by /u^ ( Frez- 



zotti et al. . 2001 1 and a spectral gap of size /z around 
the real axis (Gattringer and Solbrig 2005). Further- 



more, the twisted mass fermion matrix fulfills a general- 
ized form of the 75-hermiticity condition of the Wilson 



operator (72) 



M 



75nMi„Ti75 



(79) 



that similarly implies the appearance of eigenmodes in 
complex conjugate pairs. In contrast to Wilson fermions 
however there are no real eigenmodes due to the spec- 



(77) tral gap (Gattringer and Solbrig 2005) and therefore no 



exceptional configurations. Note that (79) also holds for 
the non-degenerate Mtm from (78). 



Numerical evidence has been found with twisted mass 
fermions for a line of first order phase transition in the 
bare quark mass that extends into the twisted mass di- 
rection ( Farchioni et al. 2005a|b I . This observation can 
be understood in terms of the phase structure of lattice 



QCD with Wilson fermions as proposed by (Sharpe and 



Singleton 1998 ) from analysis of the effective chiral po- 



tential. It represents one of two possibilities of finding a 
minimum - the other one being the appearance of an un- 
physical phase where parity and flavor are spontaneously 
broken (lAokil 119841 lAoki et al.\ 119971) . Evidence for the 



Aoki phase was found at coarser lattice spacings (Ilgen 



fritz et aL[ |2004| [Sternbeck et aZ.[|2004[ ). 

In general, it is mandatory for all simulations with 
Wilson-type fermions to avoid being too close to the line 
of first order phase transition. The situation is partic- 
ularly challenging to twisted mass fermions at maximal 
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twist however since the hne of first order phase transi- 
tion opens up in the twisted mass direction at the crit- 
ical bare mass. The minimum pion mass one can reach 
with twisted mass fermions at a given lattice spacing 
is estimated by ( |Shindler[ |2008[ ) to be ~ 300 MeV at 
a ^ 0.07 — 0.1 fm which is roughly consistent with recent 



numerical results (Baron et al. 2010a) 



For further details about twisted mass fermions we re- 



fer the interested reader to a recent review ( Shindler 



2008) 



ultralocal, i.e. they can not be realized by couplings to a 
finite number of nearest neighbors. It is therefore neces- 
sary to prove the locality of chiral fermion actions in the 
sense that the coupling decreases exponentially with dis- 
tance with an exponent that is on the order of the cutoff 
and not a physical mass. 

Different fermion operators fulfilling the Ginsparg- 
Wilson relation have been suggested. The overlap op- 
erator (Narayanan and Neuberger 1993a|b 1994 1995 



Neuberger 1998b|c I is an explicit construction that cor- 



4. Chirally symmetric fermions 

Although the Nielsen-Ninomiya theorem does not al- 
low one to retain the continuum form of chiral symme- 
try for local, doubler free fermions, an exact symmetry 
may be found at finite lattice spacing that goes over 
into the continuum chiral symmetry upon removal of 



responds to the unitary part of a Wilson operator at neg- 
ative bare mass —p. In the massless case, the fermion 
matrix is given by 



Mo 



Mw{-p) 



^mI,{~p)Mw{~p), 



the cutoff (Ginsparg and Wilson 1982 Hasenfratz et al. 



1998 Luscher 1998 Narayanan and Neuberger 19951. 



and a real mass term may be added as 



Mo{m) == Mo + Im 



(85) 



(86) 



The continuum form of chiral symmetry implies that the 
massless fermion operator M does anticommute with 75. 
One can generalize that relation by introducing a modi- 
fied 



The overlap operator ( 85 ) obeys the Ginsparg- Wilson re- 

^. Locality of the over- 



lation with an ultralocal R 

lap operator has been established numerically (iHernan 



75 = 75 (1 - 2aRM) 
and demanding that 

75 M + Af75 = 



(80) 



(81) 



dez et al. 1999) (see also (Golterman and Shamir 2003)). 
Note that by construction the overlap operator is normal. 
Overlap fermions are numerically extremely demand- 
ing. In contrast to all previously discussed fermion dis- 
cretizations, the fermion operator Mo is not sparse. The 
multiplication of Mo on a vector has to proceed via 



approximating the inverse matrix square- root in ( 85 1 . 



The condition (81) is known as the Ginsparg-Wilson rela- While it is possible to do this with polynomial (Giusti 



tion ( Ginsparg and Wilson 1982 ) and the operator R in 



(81 1 has to be local and is known as the Ginsparg-Wilson 



et al. 2003) or rational (Edwards et al. 1999 van den 



Eshof et al. 



2002 Neuberger 1998a) approximations, it 



kernel. From (80) one can see that the action 



is invariant under the chiral symmetry 



* 
* 



* (1 + ^75) 
(l + «75)* 



(82) 



(83) 



typically requires at least O(IOO) apphcations of the ker- 
nel Wilson operator to perform one matrix-vector multi- 
plication with Mo. 

Due to the exact chiral symmetry, overlap fermions 
are free of exceptional configurations at finite mass. 
This leads however to a nontrivial technical problem for 
dynamical overlap fermions that was first observed by 



(Fodor et al. 2004). The nonanalyticity of (85) implies 



The continuum form of chiral symmetry can be regained 
in observables by constructing them with the chirally ro- 
tated fermion field 



a divergence of the fermionic force term (391) at certain 
points in configuration space. Specifically, such a diver- 
gence occurs at topological sector boundaries where the 
number of zero modes changes. These points have to be 



* = 1* 



1 = 1- aRM 



(84) 



instead of the bare VP. As the full chiral symmetry is 
preserved by Ginsparg-Wilson fermions, all consequences 
of this symmetry such as the appearance of exact zero 



treated separately in the HMC integration ( Gundy et al. 
2009[ |Egri[ |2006| [Fodor et al] |2004[ ) specifically at fine 



lattices, as the simulation can get stuck in one topolog- 
ical sectorrl Alternatively, one can artificially constrain 



modes, an exactly conserved axial current (Kikukawa and 



Yamada 1999 ) and the anomalous breaking of the fiavor 
singlet part of the symmetry are also retained. It has 



been proven by (Bietenholz 



1999 Horvath 



1998) that 



chirally symmetric lattice fermion operators can not be 



® Note that while the nonanalyticity problem is particularly ap- 
parent for overlap quarks its origin is physical. Upon removing 
the cutoff, every fermion formulation should ultimately develop 
nonanalyticities at the topological sector boundaries. 
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the simulation to a single topological sector via the ad- 
dition of extra Wilson fermions with large negative mass 



( Fukaya et al. 2006 Izubuchi and Dawson 2002 Vranas 



2006 ) and treat this constraint as an additional finite vol- 



ume effect. 

Historically, overlap fermions were first formulated in 
five dimensions based on the realization that one may 
have chiral domain wall defects in a 2n -|- 1 dimen- 



From a numerical perspective, fixed point actions are 
expensive to simulate. In principle, fixed point fermion 
operators are not sparse matrices and neither can they be 
explicitly constructed out of a sparse matrix as is the case 
for overlap fermions. Consequently, one needs to truncate 
the operator to a finite range and the chiral symmetry 
is only approximate. The resulting additive mass renor- 
malization is small however and no problems with excep- 



sional vector-like gauge theory ( Callan and Harvey 1985 tional configurations have been seen (Gattringer et al. 



Frolov and Slavnov 1993| Kaplan[ 1992 ) . This five di- 2004 ) . The same paper also reports that the numerical 



mensional form of chiral fermions was further developed 



by ( Shamir 1993 ) and is known as domain wall fermions. 



Domain wall fermions do have an exact chiral symmetry 
only in the limit that the fifth dimension is large. On 
the lattice, this can of course not be realized and there is 



a remnant breaking of chiral symmetry ( Blum and Soni 



1997) (for a recent update on the size of this effect see 



e.g. (Aoki et al. 2011)) 



On a technical level, domain wall fermions are real- 
ized by 5-dimensional Wilson fermions at a negative bare 
mass M. The gauge field is 4-dimensional only and iden- 
tical for each slice in the fifth dimension. Along the fifth 
dimension s, gauge links are set to U5{x, s) = 1 for s 7^ 
except at the defect location s = 0, where they are set 
to (75(2;, 0) = Im. According to ( Shamir[ 1993), a left 
resp. a right handed chiral mode will form at the posi- 
tive resp. negative side of the defect in the limit of an 
infinite fifth dimension and these modes couple with the 
mass term m. A remnant coupling of the chiral modes 
through the bulk will appear for a finite fifth dimension 
that will be suppressed exponentially in the size of the 
fifth dimension N^. 

Due to the residual chiral symmetry breaking, domain 
wall fermions do suffer a small additive mass renormal- 
ization known as residual mass rrires- Like in the case 
of Wilson fermions it is therefore necessary in principle 
to take negative bare mass values for reaching arbitrarily 
small but positive renormalized quark masses. One could 
therefore encounter exceptional configurations, but due 
to the smallness of TOros this is not a problem in current 
simulations. The extent of the fifth dimension is typically 
around N^ = 16 in present day calculations rendering 
domain wall fermions numerically more expensive than 
Wilson fermions by about an order of magnitude. 

Another variant of chirally symmetric lattice fermions 
operators is known as perfect action or fixed point 



fermions (Bietenholz and Wiese 1996 DeGrand et al. 



1995 Hasenfratz and Niedermayer 1994). Perfect ac 



tions are obtained by following the renormalization group 
fiow of a blocking transformation to the renormalized tra- 
jectory that ends in a fixed point. Therefore, their form 
is not explicitly given but needs to be determined by fol- 
lowing the renormalization group flow. Up to truncation 
errors, the action so obtained is classically perfect in the 
sense that it has no remaining cutoff effects in the classi- 



cost is increases between one and two orders of magnitude 
compared to Wilson fermions. For a review of truncated 



perfect action fermios see ( Bietenholz 2008 ) 



Yet another variant of approximately chiral lattice 
fermions is obtained by inserting a truncated expansion 
of a general fermion operator into the Ginsparg- Wilson 



relation (81) and explicitly solving for the expansion co- 
efficients (Gattringer 2001). Numerical properties and 



cost of this variant of approximately chiral fermions are 
roughly comparable to those of the truncated perfect ac- 



tion as demonstrated in ( Gattringer et al. 2004 ) 



E. Constructing efficient regularizations 



As mentioned in sect. |II.B[ lattice discretizations of 
continuum actions are not unique. The essential step 
in discretizing a continuum action is the replacement 
of derivative terms by lattice finite difference opera- 
tors. Disregarding quantum effects, it is easy to see how 
the discretization of derivative operators can be system- 
atically improved by adding finite difference operators 
with increasing distances. For the simple case of one- 
dimensional symmetric difference operators 



Ai/(a;) 



A2/(x) 



_ fjx + a) - f{x~a) 
~ 2a 

^nx) + ^f"{x) + 0{a') 

_ f{x + 2a) - f{x~2a) 

^r{x) + ^rix)+o{a^) 



(87) 



we find discretization errors of 0{a^). In the linear com- 
bination 



A* 



4Ai - A2 



cal theory (see sect. II. E. 2 for a more detailed discussion). 



however the 0{a^) terms cancel and therefore 

A7(x) - fix) + Oia") (89) 

Generally speaking, one can systematically improve dis- 
cretized continuum operators by taking liner combina- 
tions of lattice operators and imposing conditions on the 
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C, 



C2 



FIG. 4 The four possible forms of closed gauge loops on a 
hypercubic lattice with 6 or less gauge links corresponding to 
the members of the sets Co, ■ ■ ■ ,Cz as described in the text. 



coefficients such that the continuum limit is correct and 
leading order discretization effects are cancelled. 

Finding the proper coefficients in the linear combina- 
tion is of course not always as trivial as in the example 
above. Specifically, one needs to keep in mind that when 
computing observables in a quantum theory, the clas- 
sically computed coefficients can receive radiative cor- 
rections. One can then look at a certain set of observ- 
ables and try to cancel higher order effects in them on a 
quantum level. Such a strategy has first been suggested 
by (Symanzik 1983a|b I who realized that a lattice La- 
grangian in general is equivalent order by order in a and 



g to a continuum local effective Lagrangian. 

Other strategies of finding improved discretizations are 
based on the mean field approximation or the renormal- 
ization group as will be detailed below. One common 
feature that all of the methods share is that they use the 
freedom in defining a lattice regularization of a contin- 
uum operator in order to suppress unphysical UV fluctu- 
ations in the lattice action. There is no a priori guide for 
determining which specific improvement out of the rather 
large number of possibilities is optimal. It is therefore es- 
sential to consider the potential benefits of a specific im- 
provement in relation to its computational cost. In the 
end, the optimal action will be the one that provides the 
smallest error (including all systematics) on the physical 
observables (i.e. the hadron masses) for a given amount 
of available computer time. To find a good improvement 
strategy one therefore needs to find the right balance of 
different improvements such that the overall error is min- 
imized. 



that are of 0{a'^) vanish classically. The desired contin- 
uum action has the form 



Oo = Tr(^^,F^,) 



(90) 



All together, there are three different terms that repre- 
sent possible 0{a?') corrections to this form 



Oi = Tr (D^Ff^^D^F^,) 
02^Ti-{D,F^,D,F^,) 
03^Ti-{D^F^,,D,F„,) 



(91) 



where the covariant derivative D^^ is given by (51). In 
order to reach classical improvement, we have to demand 
that these additional terms vanish. 

On the lattice, a gauge action can generically be writ- 
ten as a sum of terms of the form 



5, = /3 ^ ( 1 - iReTr(C/(7') 



(92) 



where the U{V) are path ordered products of gauge links 
along a closed path V. For the simple Wilson gauge ac- 



tion (20), which we will refer to as 5*0 in this context, 



the set of closed gauge loops Cq consists of all elementary 
plaquettes. Going one order higher, we see that there are 
also three possible forms of 6-link loops that are the mini- 
mal extensions of the elementary plaquette. The first one 
of them is the planar 2x1 loop while the other ones ex- 
tend into the elementary hypercube (see fig. |4]). Not dis- 
tinguishing between loops that differ only by rotation, we 
label the sets of loops of these different forms Ci,C2,C3. 



The expansion of the corresponding gauge actions ( 92 1 in 



terms of continuum operators (91) up to next-to- leading 
order reads (Luscher and Weisz 1985b|c ) 



-Oo + ^^0, + .. 



Si 



4 
-20o 

-20o 

'40o 



1-' 



6 



6 ^ 



(93) 



Oi 



1 



:03 



where volume sums are implied. For a general linear 
combination 



1. Gauge field improvement 



The simple Wilson gauge action (20 1 only contains 



the elementary plaquette and has discretization errors of 
Oia"^). As demonstrated by ( Luscher and Weisz 1985b|c 



Weisz 1983), one can improve the scaling by taking 



proper linear combinations of the elementary Wilson pla- 
quette and more extended gauge loops. The coefficients 
must of course be chosen such that in the continuum limit 



one still obtains the correct continuum gauge action (21 ). 
In addition however, one can demand that all corrections 



O — / ^ CiJi 



i=0 



(94) 



the correct continuum limit is therefore imposed by the 
condition 



Co + Sci + 8c2 + 16c3 = 1 



(95) 



The cancellation of all the higher order operators (91) 
may be achieved by imposing 



Co + 20ci = 



C2 



C3 = 



(96) 
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which leads to a choice of coefficients Cq = 5/3 and 
ci = —1/12. The resulting action is known as tree level 
Liischer-Weisz action. It realizes the tree level (i.e. clas- 
sical) 0{a^) improvement of the Wilson plaquette gauge 
action (20). Even before the work of Liischer and Weisz 



these coefficients have been found by ( Curd et al. 1983 ) 
with a different matching procedure. 

As mentioned above, in a quantum theory radiative 
corrections can in general reintroduce 0{a^) terms into 
observables. Generically, these corrections are propor- 
tional to g^ such that the tree level Liischer-Weisz ac- 
tion is correct up to 0{g'^a^) terms as opposed to 0{a'^) 
in the classical theory. In order to correct for these 
radiative effects, one may look at different on-shell ob- 
servables. Looking at the scattering amplitudes of mas- 
sive gluons one can construct an action with 0{g'^a^) 
scaling by modifying Cq, ci and C2 with 0{g^) terms 
(Luscher and Weisz 1985a). The resulting coefficients 
read cq = 5/3 + 0.237.g^ ci = -1/12 - 0.02521^2 and 
C2 = —O.OOAAlg^. Note that the one loop coefficients are 
explicitly scale dependent via g"^. 

It is of course possible to go beyond perturbative 
Symanzik improvement and concentrate on specific dis- 
cretization terms that are found to be important in non- 
perturbative calculations. One case of particular rele- 



vance was found by ( Lepage and Mackenzie 1993 Parisi 



1980 ) . As the gauge field is represented on the lattice in 



exponentiated form 

== 1 + igaA^{x) 



■,2„2 



'-^Ali.) + ^ 



(97) 



the theory contains vertices with an arbitrary number of 
gluons that give rise to tadpole diagrams. Although these 
are formally suppressed by powers of a, UV divergences 
in the tadpole loops apear and these terms are in fact 
scaling as powers of g^ instead. 

This problem may be adressed in mean field theory by 
redefining the relation between the lattice C/^ and the 
continuum A^. Formally writing the gauge field U^ as a 
product of an IR and an UV part 



u,=ur^ui^ 



(98) 



and replacing the UV part by its mean field value uq, one 
can now identify the physically relevant IR part C/^ of 
the lattice gauge field with the continuum field obtaining 
a relation 



Ua 



uoe 



igaAi_,{x) 



(99) 



between the lattice gauge field C/^ and the continuum A^. 
One can therefore implement tadpole improvement by 
replacing all gauge links U^ in lattice operators by [/^/uq- 
Because UV fluctuations are dominant, one can obtain a 
good estimate of uq by simply taking the expectation 



value of the trace of a gauge link in a fixed gauge or, 
alternatively, defining 



uo - [ -ReTr([/^,) 



1/4 



(100) 



In both cases uq can be either determined in perturbation 
theory or by measuring it directly as part of a nonper- 
turbative calculation. In the latter case, care has to be 
taken to determine Uq self-consistently because it is both 
a parameter in the action and an observable. 

There is a third, independent strategy of improving 
the gauge action that is based on the renormalization 
group. As already mentioned in sect. |II.D.4[ one may 
try to follow the renormalization group flow of a block- 
ing transformation in the space of all possible actions 
towards the renormalizcd trajectory. The renormalized 
trajectory is the trajectory of the renormalization group 
flow that starts from a fixed point at the critical surface 
where the correlation length diverges. Any point along 
the renormalized trajectory therefore corresponds to the 
action at a certain finite correlation length that has van- 
ishing irrelevant operator contributions and therefore re- 
produces continuum physics without cutoff effects. Ac- 
tions along the renormalized trajectory are thus called 
perfect actions. 

The exact position of the renormalized trajectory is 
elusive since the renormalization flow is not known ana- 
lytically. Strategies of finding approximations are based 
on the fact that the renormalized trajectory is attractive 
under blocking transformations as they reduce irrelevant 
operators. One should also keep in mind that a perfect 
action generally lives in an infinite dimensional space of 
couplings that has to be truncated for practical purposes 
and that a thus truncated perfect action is not guaran- 
teed to exhibit the smallest scaling violations possible 
among actions living in that subspace. 

Studying the repeated application of a blocking trans- 
formation in a truncated subspace of gauge loops. 



(Iwasaki 1983) suggested an action of the Liischer-Weisz 
form (94) but with a set of coefficients cq = 3.648, 
ci — —0.331, and C2 — c^ — 0. Also within the same 
truncation scheme, the doubly blocked Wilson (DBW2) 
action ( de Forcrand et al. 2000 Takaishi 1996 ) has been 



obtained by double blocking from Wilson configurations. 
It has the coefficients cq = 12.2704, a = -1.4088, and 
C2 — C3 = 0. Note that these coefficients do not have an 
explicit scale dependence and therefore can only cancel 
quantum effects at the scale where they are computed. 



A different strategy has been followed by (DeGrand 



et al. 1995 Hasenfratz and Niedermayer 1994 ) who ob- 



tained a classically perfect action by a saddle point inte- 
gration around g = followed by a truncation to a rather 
large set of couplings. 



18 



2. Fermion field improvement 



In the case of lattice gauge actions, the Wilson plaque- 



tte action ( 20 ) provided a unique starting point for all 



improvement efforts. In the case of fermion actions the 
range of unimproved actions is quite diverse and so are 
their major discretization effects. There are some com- 
mon improvements that positively affect all fermion ac- 
tions, but the improvement strategies are sufficiently dis- 
tinct for different discretizations so we will discuss them 
separately. We start by first discussing the improvement 
of Wilson fermions. 



The Wilson fermion action ( 68 ) has leading discretiza- 



tion effects of 0{a). These can be cancelled classically by 



adding a two- hop term to the action (Hamber and Wu 



19831 



'HW 



^h^,i2D,,^D, 



-(20-0)+ ml* (101) 



with n = X^u ^M ^'id the covariant two-hop operators 



from ( 65 ) . One could in principle compute the coeffi- 



cients of the one-hop and two-hop terms in perturbation 
theory or nonperturbatively to achieve further improve- 
ment. This was not further pursued in practice however 



since ( Sheikholeslami and Wohlert 1985 ) discovered that 



one can remove 0{a) discretization terms with a more 
local operator. This additional term is the discretized 
magnetic moment operator and the corresponding action 
reads 



S< 



sw 



Jy 



rcsw 



E*- 



F ^ 



(102) 



IJ,<v 



where the field strength F^j^{x) is usually obtained by 
taking an average of the imaginary parts of all the pla- 
quette around the point x. Due to the arrangement of 
the four plaquettes involved in the average, the additional 
term and the resulting action are commonly referred to 
as the clover term and the clover action. Tree level im- 
provement is achieved by setting the coefficient csvv = 1- 
The resulting action has discretization effects of 0{g^a) 
and 0{a^). Numerically, the 0{g^a) and the 0{a^) cor- 
rections are competing and it is not possible to say a 
priori which of these effects are dominant for a specific 
observable at a specific lattice spacing. Further improve- 



ment is possible by either perturbation theory (Luscher 



and Weisz 1996 Wohlert 19871, mean field tadpole im 
provement or nonperturbative methods |Luscher et al. 



1997). It turns out that the clov er ac tion (102) is equiv- 
alent to the Hamber- Wu action ( 101 1 up to 0{a'^) terms 



Martinelh et al. 1991 1 



in terms of rotated fermion fields ( Heatlie et al. [1991 



Although suggestions exist in the literature for im- 
proving the Wilson operator further by adding more ex- 



tended terms (Alford et al. 1997 DeGrand 



1998 Durr 



and Koutsou 2011), it is not much pursued due to the 



computational overhead. 



(l-a) 



+ 



a 
2 



FIG. 5 The principle of APE smearing displayed in the two 
dimensional case. The "thin" gauge link is replaced by a 
weighted average over the gauge link and the staples, which 
is then usually backprojected onto the gauge group. 



In contrast to Wilson fermions the staggered action 
has leading discretization effects of 0{a'^). These can be 
eliminated on a classical level by replacing the one-hop 
derivative of the staggered action ([57]) with a suitable 
combination of one-hop and three-hop derivative oper- 
ators such that 0{a^) terms cancel (Naik 1989). The 
three-hop term is used rather than the two-hop term in 
order not to interfere with the staggered flavor structure. 

The main concern for staggered fermions however is 
not the Symanzik improvement but rather the minimiza- 
tion of taste breaking effects. Since the different fermion 
components sit at the edges of the Brillouin zone, they 
interact via the exchange of hard gluons with momenta 
on the order of the cutoff scale. Suppressing these inter- 
actions, i.e. reducing the unphysical UV fluctuations, is 
therefore especially important for staggered fermions. 

The primary method in use today for reducing unphys- 
ical UV noise is link smearing (also known as UV filtering 
or fattening) . Since link smearing is used for Wilson-type 
and chirally symmetric fermion actions as well, we will 
discuss it in a more general context. Although on a tech- 
nical level it can be implemented by modifying the gauge 
fields only, it is important to remember that it strictly is a 
modification of the fermion action only. The original sug- 



gestion, put forward by the APE collaboration ( Albanese 



et al. 19871, is commonly known as APE smearing. The 



basic idea is that one can use the freedom in defining 



the parallel transport in the covariant one-hop term ( 15 ) 
of any fermion operator to suppress UV fluctuations. It 
is not necessary that one takes the same gauge link than 
used in the gauge action but instead a linear combination 
of various paths that have the correct starting and end 
points. In the case of APE smearing these paths are, in 
addition to the original gauge link, all three link connec- 
tions of the same two points (see fig. |5]) usually referred 
to as staplesjj 

One can define an APE smeared gauge link UJ}, (x) 
from the original gauge links U^{x) in d dimensions via 



t/(^P^)(a;) = (l-a)t/^(x) + 



2(<i-l 



-n^{x) (103) 



"^ Note that in the literature the sum over all three-link connections 
is also sometimes referred to as the staple. 



19 



where we used the staple sum 



(104) 



with the identity [/_^(a;) = UUx ~ fl). The smearing 
parameter a determines the relative weight of the staple 
vs. the original link and is typically set to a value a ^ 
0.6. The resulting gauge links f//i (x) are no more an 
element of the gauge group and it is therefore customary 
to backproject them onto the gauge group via 



U' 



[/(APE) 



^J7(APE)t[/(APE) 



u 



U' 



(detC/')'/^ 



(105) 



Where U' is unitary and U also has unit determinant]^ 
As the backprojection ( 105 ) is not analytic, the 



fermionic force term ( 39 1 in the pseudofermion field in- 
tegration may exhibit singularities. Although there are 
suggestions to remedy this situation by leaving out the 



second step of the backprojection ( 105 1 and use U' only 



(Hasenfratz et al. 20071, it is customary in dynamical 



simulations to use the analytic link smearing suggested 
by (Morningstar and Peardon 2004 1. They define the 
so-called stout link aqj 



V^{x) = e'''-^-~^U^[x) 



where 



with 



S^,{x) 



A^{x) - ^TtA^{x) 



M^) 



n^{x)ul{x) 



u^{x)nl{x) 



(106) 



(107) 



(108) 



The parameter p is a smearing parameter that, similarly 
to a in the case of APE smearing, determines the relative 
weights of the original link and the staple. For small 
smearing parameters, APE and stout link smearing are 
equivalent if one sets a = 2{d—l)p (Capitani et al. 20061. 
With a proper matching of the smearing parameters one 
can find a close correspondence even if their values are 



large (Hasenfratz et al. 20071. 



Both APE and stout smearing, as in general all link 
smearing techniques, generate a smeared gauge field from 
the original one, wich is usually called thin link. It is 
therefore straightforward to apply the smearing prescrip- 
tion repeatedly on the already smeared links and use this 
multiply smeared gauge link field for constructing the 



* For an alternative suggestion on doing the backprojection see 
IjDurr and Koutsou]|201l[ |. 

^ We use the term stout link as it is commonly done in the liter- 
ature. Note however that in the original paper the term stout 
link was used in a slightly different way. 
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FIG. 6 Locality of a fermion operator coupling to 6-times 
stout smeared gauge field from ( jDurr et'oL [2008 1. The stout 
smearing parameter is set to p = 0.11. At small distances 
the coupling decreases exponentially with an effective mass 
of ~ 2.2a~^ that is proportional to the lattice cutoff. At 



euclidean distances larger than v50a couplings are zero due 
to the ultralocality of the smearing procedure. 



fermionic operator. As long as the smearing parameter 
and the number of smearing steps is held constant, it 
amounts to an ultralocal redefinition of the fermion op- 
erator and does not affect the continuum limit. In fact, 
the locality range of an ultralocal fermion operator it- 
self is not at all affected by smearing the gauge links. 
Gauge link smearing does not introduce any new cou- 
plings into the fermion operator. What is affected by 
gauge link smearing is the fermion to gauge field cou- 
pling which becomes more extended. For smeared gauge 
links the fermion matrix elements are affected by changes 
of the original gauge field at further distances. If one 
keeps the number of smearing steps constant when go- 
ing to the continuum limit however, this redefinition is 
still ultralocal. In addition if the smearing parameter 
is not excessive, one expects an exponential decrease of 
the gauge field to fermion coupling within the ultralocal- 
ity range with an exponent that is proportional to the 
cutoff. This has been numerically demonstrated for a 
6-times stout link smeared tree level improved Wilson 
operator in (Durr et al. 2008) (see fig. [6]). 

A variant of APE link smearing where one tries to max- 
imize the smearing while only taking into account links 
that have a distance of at most one single lattice unit to 
the original link is known as hypercubic (HYP) smear- 
ing (Hasenfratz and Knechtli 20011. One step of HYP 



smearing combines three steps of APE smearing where 



in the first two APE steps all dimensions that would give 
distance two contributions in any one direction are disre- 
garded in forming the staple sum. The analytic version of 



HYP smearing along the lines of ( Morningstar and Pear- 
[don[ |2004j) is known as h ypercubic exponential (HEX) 
smearing ( Capitani et aZ.[ |2006 



In the context of Wilson-type fermions smearing has 
been found to drastically reduce the additive mass renor- 
malization, especially in combination with 0{a) improve- 
ment (Capitani et al. 20061. Furthermore, renormaliza- 



tion constants and the value of the improvement coeffi- 



cient csw from (102) are much closer to their tree level 



values and the normality of the operator is improved 



(Durr et al. 2005 2009 Hoffmann et al. 



2007 Hors- 



ley et al. 2008). Most importantly, gauge link smearing 



reduces the fluctuations in the real modes that lead to ex- 



ceptional configurations (DeGrand et al. 1999 Stephen- 
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staples used in APE smearing are referred to as fat3 sta- 
ples. Adding to these staples ones that extend in a second 
lattice direction and consist of five links one arrives at the 
so-called fat5 links and finally adding staples that extend 
in all three other lattice directions and consist of seven 
links each one obtains smeared gauge links referred to as 
fat 7. Finally, staples that consist of five links and extend 
two lattice spacings in one direction are referred to as 
the Lepage term ( Lepage[ [1999). Adding a Naik term 
to the staggered action, replacing the gauge links with 
fat7 gauge links plus a Lepage term and tadpole improv- 
ing the action one obtains the so-called "asqtad" action 
that is frequently used in current staggered fermion cal- 
culations and has scaling corrections of 0{g^a^). More 
recently, two steps of fat7 link smearing, the first one 
without the second one with Lepage term, with a gauge 
group projection after the first smearing step were sug- 



son et al. 2000). Thus it is possible with smeared link gested by (FoUana et al. 2007). Staggered fermions with 



Wilson-type fermions to reach physical quark masses 



(Durr et al. 2011b), which has proven to be difficult with- 



out link smearing. These are clear indications that gauge 
link smearing efficiently suppresses UV fluctuations and 
ameliorates the chiral symmetry breaking that is inherent 
in Wilson type fermions. 

One can combine nonperturbative improvement of the 
Wilson action with smearing as is done e.g. in the 
SLING (stout link improved nonperturbative clover) ac- 
tion (iHorsley et aL 1 120081). 



this variant of link smearing are known as "highly im- 
proved staggered quarks" or HISQ. 

Besides these smearing terms that were specifically de- 
signed to reduce taste splitting, simple stout link smear- 
ing has also been applied to the staggered fermion op- 
erator. Depending on the quark mass, the level of taste 
splitting has found to be generally comparable to that of 
the HISQ action or slightly less than for the specially de- 



It is also possible to use different gauge link definitions 
for different parts of the fermion operator. An opera- 
tor of this form, the so-called fat link irrelevant clover 



2010) 



signed asqtad action (Aoki et al. 2009b Borsanyi et al. 



(FLIG) fermions suggested by (|Zanotti et al.\ |2002t use ''■ Anisotropic discretizations 



an 0{a) improved Wilson operator with smeared links in 
the continuum irrelevant terms and original thin links for 
the rest. 

The improvements of Wilson type fermions carry over 
to their use as kernel operators for the overlap construc- 
tion. Overlap operators with smeared Wilson kernels are 
generally cheaper computationally, have renormalization 
constants closer to their tree level values and require less 



fine tuning of the negative mass parameter p ( DeGrand 



Kovacs 2003) 



eFari [20031 [Durr and Hoelblingj [2005al |Durr et aL[|2005| 



In the context of staggered fermions the main advan- 
tage of smeared links is the suppression of taste viola- 



tion (Blum et al. 1997 Lagae and Sinclair, 1998; Lep- 
agel |1999| |1998[ |Orginos and Toussaint[ il998, .Orginos, 



As we will discuss in sect. |III[ information about ex- 
cited states can be extracted from correlation functions 
at short euclidean distances. It is therefore desirable for 
excited hadron spectroscopy to have a fine resolution in 
time direction while keeping the resolution in space di- 
rection coarser in order to keep the overall computational 
effort small. 

Ghoosing an anisotropic discretization explicitly 
breaks the hypercubic lattice remnant of the Lorentz 
symmetry down to the subgroup of spatial cubic rota- 
tions. Gonsequently spatial and temporal lattice ex- 
tents receive different normalization and the renormal- 
ized anisotropy generally differs from the input bare one. 

Generalizing the simple case of the Wilson plaquette 



et al. . 1999|^° In the context of staggered fermions, the action ( 20 ) to include an anisotropy one obtains 



Saw = Ko E (l - ^ReTr([/,,(a;))') + ^ JI (^ " ^ReTr(C/,o(a;))] 



(109) 



For a study of reduced taste violation with an improved, un- smeared action see ^ernard et al.\^MS^ . 
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where ^o = o-o/'^o i^ the bare anisotropy factor for the gauge action, i.e. the ratio of bare spatial to temporal lattice 
spacings. Similarly one can define an anisotropic Wilson fermion operator as a generalization of the isotropic case 



(701 



M, 




(110) 



with C^ defined in ( 69 ) and v^^ the speed of light in spa- 
tial or temporal direction. In order to obtain the same 
renormalized aspect ratio for both fermion and gauge ac- 
tions one needs to tune v^ and vt. An additional tuning of 
the gauge action anisotropy is only required if one wants 
to tune to a specific renormalized aspect ratio. 

The improvement programme for gauge and fermion 
fields as outlined in sect. III. El can be carried over to 
anisotropic lattices if one keeps in mind the separation 
between spatial and temporal split and ultimate tuning of 
both gauge and fermion field anisotropies. Since spatial 
and temporal directions are in any case treated differ- 
ently for anisotropic lattices, one may even use different 
improvements on spatial and temporal field components 
that are specifically suited for either coarse or fine lat- 
tices. For further details on implementation and param- 
eter tuning of anisotropic scalar, gauge and fermion ac- 
tions see e.g. (Alford ei a\. 2001 Burgers et al\ |1988 



matescu 



et al\ 2003 ) 



Engels et al. 2000 Harada et al. 



CiieSl [20011 jCsikor et al\ |1998 



2001 



Karsch and Sta- 
19891 jKlassenl |l998l jMorrin et al. \ [2006 ( |Umeda| 



Edwards et al. 



2008 



III. EXTRACTION OF HADRON MASSES 

In sect, [n] we set up the framework for regularizing 
QCD on a discrete spacetime lattice. In this section we 
discuss how to extract the observables of interest, the 
hadron masses and energy levels, from lattice QCD. The 
emphasis in this section is on the technical details of ex- 
tracting hadron masses in a nonperturbative lattice QCD 
calculation. The question of how these measured hadron 
masses can be turned into physical predictions is dis- 
cussed in sect. IIVI 



We start by discussing the basic concept of extract- 



ing energy levels in lattice QCD in sect. III. A with an 
emphasis on the extraction of the ground state. The effi- 
ciency of this extraction of energy levels depends on the 
choice of source and sink operators, which is reviewed in 
sect. |III.B| Finally, we discuss the particular challenges 



involved in extracting excited states in sect. III.C 



A. Extraction of energy levels in lattice QCD 

The principle of extracting energy levels of hadrons 
from lattice QCD is relatively straightforward. Given 
a specific fermion matrix M{U) on a gauge field back- 
ground U, the Feynman propagator Sij{x,y) of the 
fermion field on this given gauge field background is 



Su{y,x) 



i^Iu'] 



V,x 



(111) 



where we have suppressed additional color and spinor in- 
dices that both M and S carry. These quark propagators 
on a fixed gauge field background are the basic building 
blocks from which hadronic observables may be built. 
Note that for every action that fulfills the 75-hermiticity 



condition ( 72 ) one can write 



Su{x,y) ^75Slj{y,x)j5 



(112) 



to exchange source and sink points of the propagator. 
For staggered fermions an equivalent relation is provided 



by the e-hermiticity ( 60 1 that implies 



Su{x,y) = eSlj{y,x)e 



(113) 



Having constructed a hadronic observable, one can sim- 
ply average it over the configurations that were produced 



using an importance sampling technique (see sect. II. C) 



to obtain the path integral expectation value (24) up to 



a statistical precision which is limited by the size of the 
ensemble of configurations. 

Let us assume that we are interested in the mass of 
a certain hadronic state \h) that we do not know how 
to construct explicitly. We choose two (not necessarily 
different) interpolating operators Oi/f that have a non- 
vanishing overlap with \h) 



'n/f\h)^o 



(114) 



and compute the expectation value of the correlation 
function 



G{t,0)^{0\Of{t)OlmO) 

= (0|e«*OH0)e-«*O|(0)|0) 



(115) 



between times and t. Inserting a complete set of eigen- 
states of the Hamiltonian H in the standard fashion we 
find 



G(t,0)=^ 



{0\Of\n){n\Ol\Q)^_E^ 



2En 



(116) 
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where En is the energy of the n'^ eigenstate of H above 
the vacuum energy. If the state \h) we are interested in 
happens to be the lowest energy state with the quantum 
numbers of Oi//, one simply needs to go to asymptotic 
times for all other states to die out in the correlation 



function (114) 



G(t,0) 



{0\Of\h){h\Ol\0) 
2Mh 



'Mht 



(117) 



and to extract the mass Mh and the product of matrix 
elements {Q\0 f\h) {h\0\\Q) from it. 

It is of course not possible to go to asymptotic times on 
a finite lattice. However, since the higher energy states 
are dying off exponentially in euclidean time with an ex- 
ponent that is their energy difference to the ground state, 
it is often possible to reach a distance that is effectively 
asymptotic for limited lattice time extents. The time 
interval that starts from where one can not see the ef- 
fect of higher energy states and ends at a possible loss 
of signal is called the plateau region. For ground state 
spectroscopy it is desirable to extend this plateau region 
in order to get a statistically clean signal. This is often 
achieved by choosing the operators Ci// in such a way 
that they have a large overlap with the ground state and 
small overlap with all other states. This will further be 
detailed in sect. IIII.BI 

Another common strategy to extend the plateau range 
is to take cither Oi or Of to be a sum of local operators 
Oi over an entire time slice 



0,,f{t) = Y,Oi{t, 



(118) 



With this choice either the initial or the final state is 
projected to zero spatial momentum and consequently 
all higher momentum excitation that may appear in the 



sum over states (116) are cancelled 



Since our lattices have a torus topology with a period 
T in time direction, there is a backward contribution that 
is dominant for T— t > t. It has a similar form than ( 117 ) 



with the difference, that the complete set of states has 
now been inserted on the other side 



G(t,0) 



T-j^oo A0\Ol\h){h\Of\0) _ 
2Ms ^ 



Mj^iT^t) 



(119) 



In addition one also has in principle contributions from 
propagators that wrap one or more times around the lat- 
tice in time direction. These contributions are tiny how- 
ever - each additional wrapping gives a suppression factor 
g-TMh Qj. ^-TMf^ _ g^j^^ |.]-^g resulting geometric series can 
be summed up. Putting all this together, we find that in 



the plateau range the correlation function (115) is given 
by 



G{t,0)^Af 
+Ab 



o-Mht 



^~Mf^(T~t) 

2Mf^{l-be-TMn) 



(120) 



with the matrix elements 



Af^{0\Of\h){h\ol\0) 
At = b{0\ol\h)(Ji\Of\0) 



(121) 



One undesirable feature of ( |120 ) is the exponential de- 
cay of the signal with euclidean time. For large time 
separations t or T — i, the signal exponentially vanishes. 

In order to check for the existence and extent of the 
plateau region, one can define an effective mass 



Moff(^ -I- a/2) = In 



G{t + a,0) 
G{t,0) 



(122) 



which will be Mh or —Mh' in the region wher e eit her 
the first or the second exponential dominate in ( |l20[ ) p^ 
As one can see in fig. [7j the effective mass plot is very 
useful for identifying the plateau region of a correlation 
function. One should however keep in mind that the 



time t for which the asymptotic regime (117) is reached 



may vary widely. It is possible that the coupling of an 
operator to the lowest energy state of the same quantum 
numbers is nonzero but so tiny that the ground state is 
not reachable. 

We now turn to the explicit form of the interpolat- 
ing operators Oi/f. For Wilson fermions the simplest 
form they can assume is that of a local operator with 
the correct quantum number. For pseudoscalar mesons 
composed of two different quark flavors one can e.g. take 



Pi(x) = ^i(a;)75^2(a;) 



(123) 



Note the appearance of the ground state h. It coincides 
with h except for cases where the lowest state that cou- 
ples to both Oi and Of is different from the lowest state 
that couples to both 0| and 0\. The factor 6 = ±1 has 
been inserted to account for a sign flip that occurs for in- 
terpolating operators with an odd number of quark fields 
when the timeslice is crossed that incorporates antiperi- 
odic boundary conditions in time direction. Without loss 
of generality, we assume here that this time slice is tra- 
versed in the backward contribution. 



Aoi{x) = ^l{x)-/0-/5'i'2ix) 



(124) 



In order to demonstrate how to construct a proper lat- 
tice observable, we take as an example Of{t) = Vi{t,x) 



^^ In the case where G{t,0) is either symmetric or antisymmetric 
one can modify | |122| l such that it gives M^ff = M/^ throughout 
the entire plateau region. Sec ( Fle ming et al.\ |2009| for a more 
thorough discussion of effective masses. 
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FIG. 7 Plot of the effective mass Moff(t + a/2) vs. t/a { 122 1 



of a decuplet baryon from a recent lattice calculation ( Aoki 



et al. 20111. One can see clearly the onset of the plateau 
and the eventual loss of signal. The plateau region is indi- 
cated together with the value of the mass obtained from a 
fit to the correlation function. Plot reproduced with friendly 
permission of the RBC-UKQCD collaboration. 



and Oi{0) — Vi{0, 0) and plug these operators into ( 115 1 
Using X — {t,x), the resulting Greens function is 



Gpp{x,Q)^{Q\Vi{x)Vl{m) 

= (0|§i(i)75*2(a;)*2(0)75*i(0)|0) 



= (0|*i(a;)75*2(a;)*2(0)75*i(0)|0) 
= -{TY{Si{Q,x)j,S2{x,0h5)) 



where in the last line the expectation value (• • ■ ) denotes 
the (properly weighted) average over all gauge configu- 
rations and Si denote the Feynnian propagator of the 
quark flavor i on a given gauge field background. A 
graphical representation of this greens function in terms 
of quark propagators is displayed in fig. (Saj Using the 
75-hermiticity relation (112), (125) can be cast in the 
form 



G 



pp 



{x,{)) ^ - (ty (s\{x,0)S2{x,Q) 



(126) 



which can be obtained in practice by computing the in- 
verse of the fermion matrix on just one source point 0. 
Note that this is true for an arbitrary sink point x so 
that in particular one does not need to perform more in- 
versions when either summing over all sink points a? in a 
given time slice or computing the Greens function from 
to a different sink point in an arbitrary time slice. 

In case of flavor singlet interpolating operators 



V{x) = *(a;)75*(x) 



(127) 




(c) Flavor singlet disconnected 





FIG. 8 Graphical representation of the contraction for a fla- 
vor non-singlet I (a) I mesonic Greens function (|125|126 1 and of 



the connected | (b) | and disconnected (c) contributions to a fla- 



(125) vor singlet mesonic Greens function (128 1 



the Greens function contains one more Wick contraction 

G%{x,Q)^{0\V{x)V\m) 



|*(a;)75*(x)§(0)75*(0)|0) 



|*(a;)75*(x)*(0)75*(0)|0) 



(128) 



+ (0|^(x)75*(a;)^(0)75*(0)|0) 
= -(Tr(5t(x,0)5(x,0))) 

+ (Tr (5(0, 0)75) Tr(5(a;,a:)75)) 

that leads to a quark line disconnected contribution to 
the Greens function that is displayed in fig.lSclin addition 



to the connected piece (fig. 8b) that is also present in 
the flavor non-singlet case ( |126 ). Since the source and 
sink points coincide in the disconnected piece, the 75- 



hermiticity relation (112) does not provide any further 



simplification as in the case of the connected piece. 

For twisted mass, chiral and clover fermions the con- 
struction of operators has to be done in the properly 
transformed basis (see ( 76|84 ) for the case of twisted mass 
and chiral fermions respectively) but is identical to the 
Wilson case otherwise. 

For staggered fermions the construction of interpolat- 
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ing operators is complicated by the presence of four inter- 
acting tastes for each fermion flavor. Roughly speaking, 
when using the spin- flavor basis introduced in sect. |II.lXT] 
as a guide, one can construct mesonic interpolating op- 
erators along the same lines as for Wilson type fermions 



(the rigorous construction can be found in (Golterman 
and Smit| |1 9 8 4b; Kil cup and Sharpe| |1987| |Klube7g 
Stern et al. 1983j)). One should note however that due to 
the distribution of spin degrees of freedom within an ele- 
mentary hypercube all operators that have a taste struc- 
ture different from their spin structure are necessarily not 
localized to a single lattice point. One relevant example 
of an operator that is local and also leads to a correlation 
function that is positive on every gauge configuration is 



v!iy)^Xiiy){j5^^5)X2iy) 



(129) 



for which the goldstone pion is the ground state. 

The simplest baryon operators for Wilson quarks that 
have the nucleon as a ground state are the local 



A/"/ = eabc (ulCj^db) Uc 



and 



^fl = eabc (ulCdb) 75Wc 



(130) 



(131) 



where a, b, c are color indices and the charge conjugation 
operator C = 7072 s. We have suppressed explicit spin in- 



dices and coordinates. The difference between ( 130 ) and 



( 131 1 is that in the first case one starts off with a pseu- 



doscalar "diquark" (the bracketed expression) whereas in 
the second case the diquark is scalar and the 75 only en- 
ters when combining the diquark with the remaining u. 
This difference in spin structure leads to a very different 
nonrelativistic limit of the two operators - 0(1) for J\f^^' 
vs. 0{p^/E^) for A/"'^^ - which in turn implies that the 
relative overlap of A/"*-^-* with the ground state as com- 
pared to_the_excited_statcs is much larger than that of 
A/'(2) ( |Bowler et aT]|T984j Leinweberj_1995} . For decuplet 
baryons, one can construct an interpolating operator by 



replacing the pseudoscalar diquark in (130) with a vec 



tor one (Chung et al. 1982, Leinweber et al. 1992). An 



interpolating operator coupling to the A++ can for ex- 
ample be obtained by 



2?p = eabc {u^Cjf_,Ub) Uc 



(132) 



Greens functions G^^ that arise from using ( 132 ) for both 



source and sink are however not pure spin 3/2 (jLeinwe- 
ber et al. 1992). A spin projection to a pure spin 3/2 
state may be performed using the projection operator 



( Benmerrouche et al. 1989[ Van Nieuwenhuizen[|198l") 



PW^ ^ -^^i^- 37^7^-7^ (7rPr7MPM +PM7i^7rPr) (133) 

Numerical evidence suggests that even the unprojected 
correlator couples almost exclusively to the spin 3/2 state 



(Leinweber et al. 1992) 



For staggered fermions, zero momentum baryon oper- 



ators have been constructed by (Golterman and Smit 



1985). This construction does not rely on the spin- flavor 



interpretation of baryons but instead is based on analyz- 
ing the discrete spacetime symmetry group of staggered 
fermions. One interesting feature of this construction is 
that there exists an operator that couples to the A but 
does not couple to the nucleon while on the other hand 
every operator that couples to the nucleon also couples to 
the A. Since the mass difference SM of these two states 
is rather small, the A dies off slowly in euclidean time 
t as ex e~**^* and it is challenging to find a plateau re- 
gion. Effects of different quark flavors have recently been 



added to this formulation by ( Bailey 2007 ) . 



An alternative approach to extracting hadron masses 
via measuring the free energy of a system with a fi- 



nite baryon density has been suggested by (Fodor et al. 



2007). This method however requires the introduction 



of a chemical potential and its practical use is severely 
limited due to the sign problem. 



B. The role of operators 



As we have seen in sect. |III.A[ the choice of operators 
determines the relative strength of the coupling to dif- 
ferent states with the same quantum numbers and there- 
fore potentially has a huge influence on the quality of 
signal. We first consider the case of the simple method 
of improving a local operator by summing it over either 
source or sink timeslice to obtain a projection to zero 
spatial momentum. 



When writing the hadron correlation function (115) in 
terms of Wick contractions of the quark fields, it is of- 



ten possible to use the relations ( 112 ) or ( 113 ) in order to 



obtain an expression where the quark field sources are re- 
stricted to timeslice (an explicit example for the flavor 



non-singlet pseudoscalar propagator was given in (126)). 
In these cases, if one only needs to do the momentum 
projection at the sink timeslice i, it is a trivial summa- 
tion. Performing the momentum zero projection on the 
source side in principle requires computing the inverse of 
the fermion matrix for every point in the source times- 
lice t = 0. Since this would involve a prohibitively large 
number of fermion matrix inversions, a range of meth- 
ods has been developed to deal with this problem and 
related ones where one needs information of the fermion 
propagator from every source point to every sink point 
on the lattice (Bali et al. 2010 2005 Bernardson et al. 



T99m 



and t 



Bitar et al. 



eller 2002 



19891 IBoucaud et al.\ |2008[ [DeGraiidl 
de Divitiis et al.\ 1996 Dong and Liu 



1994[ [Duncan and Eichten[ |2002| [Eicker e^ aL[|1996[ [Fo- 



ley et al.\ I2005[ [Kuramashi et al.\ |1993[ [McNeile and 



Michae l, 2001; Michael and Peisa[ [T9981 [Neff et a/.[[200T 



Wilcox[[1999[ ). These methods, commonly known as all- 
to-all techniques, are based on either stochastic estimates 
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or eigeniTLode approximations of the full propagator or a 
combination of both. 

The basic idea behind a stochastic estimate of the all- 
to-all propagator is to compute the inverse of the fermion 
matrix on a number of source vectors ^^^' 



.(*) 



Mr/c^'^ 



Provided that the ^'^^^ fulfill the condition 



(134) 



J2i'-'^\y)C^'\^)=S.,yl (135) 



a single lattice point. Hadrons however are extended 
objects with quark distributions that have finite widths 
and different shapes. One strategy to improve the over- 
lap between an interpolating field operator and a specific 
hadronic state therefore consists of trying to model its 
quark distribution. One possibility that is often used is 
to replace the quark point-sources "^{t^x) by a sum over 
quark sources at neighboring sites. In its simplest form, 
one can take 



vI/'(i,y) = ^G(y,f)vI'(i,x) 



(137) 



one can obtain the propagator between any two lattice 
points on a fixed gauge background U as 



%(2/,x) = 5]C«^(y)a«(a:) 



It is often useful to restrict the vectors ^^'^ to a subspace 
of the lattice, i.e. to modify the relation ( 135 1 so that it 



is zero outside a certain subspace. If one wants to obtain 
full all-to-all propagators it is then of course necessary to 
have more than one set of source vectors ^^*^ so that all 
sets together cover the entire space desired. This method 
is sometimes referred to as partitioning or dilution and 
commonly used subspaces include individual spin com- 
ponents and individual time slices. 

We are interested in the case where th e sou rce vectors 
^(*^ are random vectors and the relation ( 135 1 is fulfilled 



stochastically. Examples of sets of stochastic source vec- 



tors that fulfill ( 135 ) are Z2 noise where every component 
of ^'•*' is randomly chosen to be either of ±1 or C/(l) noise 
where every element is a complex random number with 
unit modulus. The question of how many source vec- 
tors are necessary to get an optimal signal for a specific 
observable at minimum computational cost has to be an- 
swered numerically. It is important to note however that 
since one is usually not interested in the propagator on a 
specific gauge configuration U and the path integral sum 
does commute with the sum over all source vectors ^^*\ 
the optimal number of sources per configuration might 
even turn out to be one. 

Eigenmode approximations generally assume that the 
inverse of the fermion operator is well approximated by 
restricting it to a subspace that is spanned by a num- 
ber of its lowest eigenmodes. These eigenmodes are then 
computed and an approximation of the full matrix inverse 
can be found. This truncation represents an uncontrolled 
approximation, but it may easily be supplemented with a 
stochastic estimate of the effect of the higher eigenmodes. 
One only needs to project out the components of the orig- 
inal source vector that lie in the space of low eigenmodes, 
treat them exactly and use stochastic estimators on the 
orthogonal compliment. 

Until now we have only discussed (sums of) local oper- 
ators, i.e. operators where all the quarks originate from 



where G{y, x) is the smearing kernel that determines the 
relative weight of the source points. A computationally 
convenient restriction of the smearing kernel is the factor- 



1991) 



(136) ization ansatz (Bacilieri et al. 1988 DeGrand and Loft 



G{y,x) =g{x)g{y) 



(138) 



where effectively the quark fields are independently 
smeared. Typical choices for the form of g{x) include 
a wall (Bitar et al. 1990a), a hard sphere or box (Ba- 



cilieri ei a/. TT988 1 , a gaussian ( DeGrand 1998 DeGrand 



et al. 1998) or a radial exponential ( Ali Khan et al.] 
2OO2I) of different size. Note that in general \E''(i,y) in 



(1371 is not gauge invariant and one therefore needs to 
work in a fixed gauge on the source timeslice. A more 
sophisticated method of smearing the quark source that 
is gauge invariant is known as Wuppertal smearing (or 
gauge invariant gaussian smearing) (Gusken 1990) where 
a covariant laplacian is added to the unit operator and 
repeatedly applied N times 



G = 




(139) 



The one-hop term Vi is given by (15) and the smearing 



parameters a and N determine the form and width of 
the source that is approximately gaussian for large N 
on a trivial gauge fieldj^ In order to suppress gauge 
noise, one can use smeared links (see sect. II. E. 2) for 
the gauge field in (139). A variant of this method is 



the Laplace-Heaviside (LapH) smearing (Peardon et al. 



2009 ) where the smearing kernel is defined as a Heaviside 



step function on the lowest eigenmodes of the covariant 
Laplacian. This method requires the explicit inversion 
of the fermion matrix on a number of lowest eigenmodes 
of the covariant Laplacian that grows with the volume 
of the system. This growth with volume of the number 
of required fermion matrix inversions can be countered 



12 Sec (Allto n et aE] |l993|[Burch et a/. '20041 |Lacock et aL|[l995| l 
for different constructions of nonlocal gauge invariant operators. 
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by not computing propagators on all eigenmodes exactly 
but instead using stochastic techniques, similar to the 
ones described above, in the eigenspace of low modes to 



by (Babich et al. 2006 2007). They considered at the 



estimate them ( Morningstar et al. 2010 20111. 



Replacing the original quark sources in local hadron 
operators such as ( 123|124|130|131 ) with smeared ones 
usually does improve the overlap with the desired state. 
Especially in the case of excited states however, it is often 
useful to go further. When trying to find an operator that 
has maximal overlap with a certain state, one can use its 
quantum numbers and expected wave function to model 
a lattice operator. The two quantum numbers of interest 
are spin and parity. Parity is not broken by the lattice 
regularization and therefore one can construct operators 
that couple only to states of a given parity exclusively. 



Note however that the backward contribution in (120) 



will be of opposite parity than the forward one. For in- 
terpolating operators with an odd number of quark fields, 
it is possible to use the relative sign flip in the backward 
amplitude Ab (121 1 between periodic and antiperiodic 



boundary conditions to eliminate these. ( Csikor et al. 
20031 ILeinweber et al.\ [20051 [Sasaki and Sasaki[ |2005 



Sasaki et al. 2002 ) 



Spin on the other hand is the quantum number cor- 
responding to the continuum rotational symmetry group 
SU{2) that is broken down to the symmetry group of cu- 
bic lattice rotations. This group, the octahedral group O 



has been studied in ( Johnson 1982 ) and it was found that 



there are five irreducible representations corresponding to 
integer spin and three corresponding to half integer spin, 
all of them containing a whole tower of partially over- 
lapping continuum spin representations. One can con- 
struct lattice operators that transform irreducibly under 
the octahedral group O ( Basak et al. 2005a|b 2007 1 and 
are especially beneficial for the extraction of highly ex- 
cited baryon resonances. In addition to the quark source 
smearing, these operators generally involve quark sources 
that can each be covariantly displaced from a reference 
point by one lattice unit in any direction. The general 
form of such a baryon operator is an appropriate linear 
combination of terms of the form 



B = e,bc{D',^r{D'^-^t{D',^f 



(140) 



where ahc is color and we have suppressed spin and flavor 
indices. The operators D[ are gauge covariant displace- 
ment operators by one step in direction i (or the unit 
operator if i = 0). For more detailed reviews on the 
construction of operators for excited state baryon spec- 
troscopy see e.g. ( [Basak et al. 2006 Lang 



weber et al. 2005 1 



2008 Lein 



A different approach of finding an operator that has 
optimal overlap with the ground state was developeq^ 



See ( [Draper and McNeile[|l994| for a similar approach for heavy- 
light systems. 



sink side a general meson operator of the form 

M(t, r) = *i(i, x)T^2{t, v) X 5{\y- x\ - r) (141) 

and studied the profile in r of the resulting correlation 
function as a function of t. Once a plateau is reached, 
the profile was found to settle into an asymptotic form 
(j){r) which can then be used to construct an optimal sink 
operator 



O(i) = $^0(r)X(<, 



(142) 



A special problem occurs when one tries to extract 
the energy level of a scattering state or a bound state of 
hadrons with a non-minimal number of valence quarks. 



As repeatedly noted in the literature (see e.g. (Bulava 
et al. 2010 Engel et al. 2010)), the coupling of sin- 



gle hadron operators to multi hadron states is extremely 
small. A prominent example is the p resonance which 
has a TT — TT scattering state as a ground state in infi- 
nite volume. A closely related phenomenon occurs in 
heavy quark physics where string breaking between static 
quarks can only be observed when one introduces explic- 
itly a "broken string" operator that consists of two sepa- 



Sommer 1998) 



rated static-light mesons ( Bali et al. 2005 Knechtli and 



Both phenomena suggest that it is very difficult to pro- 
duce a pair of appropriate sea quarks from the vacuum 
that allows for the propagation of an intermediate mul- 
tihadron state. Intuitively this is understandable as the 
occurrence of large sea quark loops is suppressed in the 



path integral ( 24 ) 



For spectroscopy in channels where there is a scattering 
state below the relevant resonance or for the spectroscopy 
of exotic objects that consist of a non-minimal number 
of valence quarks - such as tetraquarks or pentaquarks - 
one therefore needs to use appropriate interpolating op- 
erators that correctly reflect the valence quark structure 
of the desired object. 

When constructing interpolating operators one usually 
tries to avoid situations where a quark line may start and 
terminate at the same time slice since these operators are 
known to be much noisier than operators where all the 
quark lines run from the source to the sink timeslice. In 
some cases however, e.g. for isoscalar mesons or generally 
for the multihadron operators mentioned above, quark 
propagators that attach to one time slice with both ends 
are unavoidable. As an example of how a combination of 
the above mentioned techniques can lead to a decent sig- 
nal in these notoriously difficult channels fig. [9| displays 
one current determination of the correlation function and 
mass plateau in the rj channel that was obtained using 
the stochastic LapH technique described above. (For sim- 
ilar results with a different mix of techniques see also 



(Alexandrou et al. 2010)) 
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FIG. 9 Correlation function (left panel) and effective mass (right panel) in the rj channel using a stochastic LapH source. 
Connected (fwd) and in the left panel disconnected (smt) contributions are plotted additionally. Plot taken from ( Morningstar 
\et al.\ |2011[) with friendly permission of K. J. Juge. 



Finally, there are also suggestions to completely avoid 
the problem of computing disconnected diagrams by 
finding relations between them and disconnected dia- 
grams in partially quenched chiral perturbation theory 



(Delia Morte and Juttner 2010). 



C. Extracting multiple energy levels 



In sect. |III.A| we have seen how to extract the ground 
state of a channel by going to asymptotic euclidean time. 
In order to extract excited state masses, one can in prin- 
ciple fit the correlation function to a multiexponential 
form 



G(t,0) = ^A 



„-A/„t 



(143) 



where we have ignored backward contributions. A fit 
of the form ( 143 1 with free parameters An and Af„ can 



typically not be stabilized numerically for more than 
two states and even extracting reliable first excited state 
masses with this technique is challenging. 

One can overcome these problem by using a variational 
method jBlossier et aL[ |2009[ [Luscher and WoI5| |1990 



Michael 1985 ) . The basic idea is to expand the basis to 
include N initial and final state operators Oi^ and Of^ 
with the same quantum numbers that ideally couple to 
different energy states preferentially and then construct 
the complete cross-correlator matrix 



GUt,o) ^ {o\Of,it)oUom 



(144) 



One can then define a matrix M(t,to) from 

G(t,0)=M(i,to)G(io,0) (145) 

and analyze its eigenvalues A„(t,io) Sind eigenvectors 
Vn{t,to)- At large euclidean times t and to the eigen- 
basis of the matrix M will then align with the eigenbasis 
of the Hamiltonian and the eigenvectors will behave as 



A. 



-B„(t-to) 



(146) 



from which one can extract an effective mass for each of 
the energy levels similar to ( 122 ). The energies of a num- 



ber of lowest lying states in a given channel can thus be 
determined provided that the operator basis chosen has 
sufhcient overlap with all of these states and the quality 
of the data is good enough do determine all the elements 
of the transfer matrix with sufficient accuracy. 

A robust variant of the variational method was sug- 
gested by ( Mahbub et al. 2009a ) ^ Instead of extracting 
the energies from the eigenvalues directly one can use the 
elements of the eigenbasis to project out effective single 
state components from the matrix valued propagator and 
analyze them with standard single channel methods. 

An example fig. [T0| displays masses of the ground state 
and first three excited states of the nucleon channel ex- 
tracted with the variational projection method. 



See i JDraper et al.\ |l995[ l for a similar method for heavy-light 
systems. 
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FIG. 10 Plot of the masses of the ground and excited states 
in the nucleon channel. The extraction has been performed 
by fitting single state correlators that were obtained by pro- 
jection of a 4 X 4 matrix correlator onto the elements of the 
eigenbasis of the transfer matrix M(t, tp) for different t and 



to- Plots are taken from (Mahbub et al. 2010c I with friendly 
permission of D. Leinweber. 



IV. PHYSICAL PREDICTIONS 

Having extracted hadron masses from nonperturbative 
lattice QCD calculations, we now need to turn them into 
physical predictions. 

QCD, and also lattice QCD, is a theory without any 
intrinsic scale - it is formulated in terms of dimcnsionless 
quantities entirely. In order to extract from it dimcnsion- 
ful quantities such as hadron masses, we need to supply 
it with a scale. We can do this in general by picking 
a scale setting observable, a quantity that is dimension- 
ful in the continuum theory. In the lattice theory we can 
then measure the appropriate dimcnsionless combination 
of any target observable we wish to extract with the scale 
setting observable and get a dimensionful prediction for 
the target observable by fixing the scale setting observ- 
able to its dimensionful input value. One such target 
observable would be the lattice spacing a itself. 

The scale setting procedure outlined above is not en- 
tirely sufficient for making physical predictions. Once 
we have supplied lattice QCD with a scale we still need 
to fix its remaining parameters - the bare quark masses. 
For Nf non-degenerate quark flavors we can generally 
do so by fixing Nf linearly independent dimcnsionless 
observables in the lattice theory to their desired input 
value. What one would ideally like to do then is to fix 
the Nf + 1 dimcnsionless bare parameters of the lattice 
theory, the bare quark masses and the gauge coupling, 
such that the Nf dimensionless observables on the lat- 
tice assume their physical values exactly and the lattice 
spacing a is of the desired size. One could then measure 
any observable on the lattice for a range of lattice spac- 



ings a and, with the appropriate functional form that is 
given by the discretization effects of the specific action 
used, extrapolate them into the continuum a = 0. 

There are a few obstacles towards implementing this 
ideal procedure in nonperturbative lattice calculations. 
The first one being that one can not simply go to the 
physical point by setting the bare coupling and input 
quark masses in a lattice calculation to their physically 
observed values since these physical values can not di- 
rectly be measured in experiment. Quantities that are 
accessible experimentally, such as hadron masses or de- 
cay widths, have relations to the bare parameters that 
need to be determined on the lattice themselves. One is 
therefore left with the choice of either tuning the bare pa- 
rameters or computing both target observables and those 
used for scale setting and parameter fixing at various un- 
physical points followed by an interpolation to the phys- 
ical point. 

When trying to implement either of the two proce- 
dures, tuning or interpolating to the physical point, one 
faces the more technical problem that it is very difficult 
to reach the physical point for light quarks. As discussed 
in sect. |II.D| there is a variety of different fermion dis- 
cretizations but each one does have a specific problem 
when making the quark mass light. In the case of Wil- 
son fermions one is faced with exceptional configurations, 
staggered and twisted mass fermions have the problem of 
flavor/taste splitting and chiral fermions are simply ex- 
pensive in terms of computer time needed. None of these 
problems is insurmountable but they turn out to be suffi- 
ciently severe so that an extrapolation in the light quark 
mass to the physical point is still the rule rather than the 
exception in hadron spectroscopy calculations. 

A third problem that is less severe in practice but still 
has to be considered is that nature is not QCD alone, 
even for the light hadron spectrum. Quarks are electri- 
cally charged and there are QED corrections to hadron 
masses. Similarly, isospin symmetry is usually assumed 
in lattice calculations but broken in nature. Both of 
these effects are relatively small as can be seen from the 
isospin splitting of the experimentally observed hadron 
spectrum, but they still need to be considered. 

We start in sect jIV.A by reviewing the problem of 
reaching the physical point. In sect. |IV.B| we discuss 
how to remove the cutoff i.e. how to reach the contin- 
uum limit for the various lattice discretizations. Finally, 
IV.C[ we consider finite volume corrections that 



in sect. 



are especially relevant for resonant states and we con- 
clude with the discussion of subleading effects from QED 
and isospin breaking in sect. |IV.D| 



A. Reaching the physical point 

Having extracted hadron masses from simulations we 
now need to turn them into physical predictions. As 



29 



discussed in the beginning of sect. IV we need to go or 
extrapolate to the physical point and we need to set the 
scale. 

On a technical level tuning to the physical or any other 
target point is achieved through tuning the bare param- 
eters of the lattice theory - the coupling f3 and the bare 
quark masses in the lattice Lagrangian. Defining the 
physical or any other target point on the other hand is 
generally done by comparing dimensionless combinations 
of continuum observables with the corresponding lattice 
observables and the scale can be set by comparing any 
one dimensionful continuum observable. Provided that 
all effects beyond QCD with the given number of flavors 
are correctly accounted for and provided also that the 
chosen lattice discretization of QCD does have the cor- 
rect continuum limit, all possible choices of finding the 
physical point are equivalent in the continuum limit. If 
any of the above assumptions is violated, such as is the 
case e.g. in quenched QCD, continuum predictions are 
not unique and depend on the specific choice of defining 
the physical point. 

It is usual to treat the two light quark flavors u and d as 
degenerate and to include isospin breaking as well as elec- 
tromagnetic effects as corrections. One can then define a 
light quark mass rh = (m„ + md)/2. In order to tune to 
the correct light quark mass one combination that is often 
used is the ratio of the pion mass, the square of which is 
proportional to the light quark mass at leading order, to 
an observable that depends less on the light quark mass. 
In early quenched work on lattice hadron spectroscopy 
the QCD string tension was often used for this purpose 



in lattice hadron spectroscopy. 

At the physical point the scale can then be set by com- 
paring a dimensionful continuum observable with the cor- 
responding lattice observable. At any other non-physical 
point the scale setting is conceptually ill defined. It is 
nonetheless usual to either use the scale defined at the 
physical point for all theories with the same coupling /3 
and arbitrary quark masses (mass independent scale set- 
ting) or to obtain the scale by comparing a dimensionful 
continuum observable with its lattice counterpart at any 
nonphysical point (mass dependent scale setting). 

In computations with a dynamical strange quark its 
mass also has to be set. In analogy to the light quark 
case, this is usually done via the kaon mass, the combi- 
nation M|- — M^/2, which is proportional to the strange 
quark mass to leading order, or the mass of the ficti- 
tious 77s , the pseudoscalar meson that is composed of two 
strange valence quarks. 

As for the specific functional form of an extrapolation 
or interpolation to the physical point the most straight- 
forward form is one that is linear in the quark mass. Since 
M^ (X rh and Mj. = M]^ - M^/2 oc rUs to leading order 



(Gell-Mann et al. 1968), one can fit any other hadron 



mass Mx to a form 



Mx=a + hM; + cM\ 



K 



(147) 



Going beyond this leading order one can perform sys- 
tematic expansions around either a point with two or 
three massless quark flavors or a general, massive point. 
For the first case, chiral perturbation theory (xPT), an 
effective field theory based upon the chiral symmetry pat- 



( Bernard ei a/. 1983 Creutz[ 1980a Creutz ei aL 1983 tern of QCD has been developed ( Gasser and Leutwyler 



et al 



nen[|1981[|Weingarten[|1982| as lat er on was Mp ( [Fucito 

In fui: 



Hamber and Parisi[ 1981[ Marinari et aL[|1981a[ Pietari 



19821 or M$ (Lipps et al. 19831. 



1984 1985 Weinberg, 1979). It provides an asymptotic 



QCD 

these choices are not optimal. The vector meson p is a 
broad resonance and not the ground state in its chan- 
nel, the singlet ip contains disconnected diagrams and 
the string tension is ill defined with dynamical quarks. 
Even in the quenched approximation any resonant state 
mass is not an ideal choice for a scale setting observable 
as the connection of the measured ground state mass to 
the experimentally measured mass above decay thresh- 
old is not obvious. Quantities that are frequently used 
today are either baryon masses that are stable in QCD 
(Mtv, Mh, Mq) ([Alexandrou et al 



expansion around either the two or three flavor massless 
point. As it is built around the assumption of spon- 
taneous chiral symmetry breaking in the massless the- 
ory, it is particularly suited for treating properties of the 
pseudo-Goldstone bosons of this symmetry, i.e. the pions 
and, to a somewhat lesser extent, the kaons. 

For masses other than the pseudoscalars, xPT generi- 
cally predicts a leading nonanalytic term of the form M^ 
(Langacker and Pagels 1974). More formally, one can 



include baryons in xPT ( Gasser et al. 1988[ ) but the re- 
sulting series is only slowly converging. An alternative 
formulation with better convergence properties is heavy 



et al.\ 2008 |Engel et al. 2010 



masses of baryon multiplets (Bietenholz et al. 2010b) 



Lin et al. 20091, average 



2009 Aoki et al. baryon xPT ( Bernard ei aL 1992 Jenkins and Manohar 



2009a[ 2010[ [Brandt et al. 2010 Bulava et al. 2010 Durr 1991) which treats baryons as nonrelativistic particles 



and currently is most commonly used to fit lattice baryon 
data. An extension of heavy baryon xPT for staggered 



or distance measures (ro, ri) in the heavy quark poten- 
tial (Sommer 1994) that are in turn determined from T 



fermions was developed by (Bailey 2008). Recently, the 



covariant approach (Becher and Leutwyler 1999) that 



spectroscopy (Bazavov et al. 2010a Bernard et al. 2001 promises better convergence behavior for heavier pion 



Davies et al. 2010). The use of matrix elements, such as masses has been revived (Dorati et al. 2008 Durr et al. 



the pseudoscalar decay constants, although common in 2010 2012 1 and used for chiral fits of the baryon octet 



other areas (Blossier et al. 2009 Borsanyi et al. 2010 



Giusti et al. 2001 Noaki et al. . 2009 1 is less often seen 



Partially quenched heavy baryon X-PT is an exten- 
sion of heavy baryon xPT where the sea and valence 
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quark masses may have different values ( Beane and Sav- 



age 120021 Che^^nTS^^gel 2002 Labrenz and Sharpe 



1996[ Savage 2002). It is useful to describe lattice data (Aoki et al. 2010), the PACS-CS collaboration has 



In order to avoid these problems, reweighting tech- 
niques have recently been applied to this problem. In 
the 



with multiple valence quark masses for each sea quark 
mass and for describing "hybrid" calculations where the 
sea quarks are regularized differently than the valence 
quarks. This technique is sometimes employed to keep 
the computational costs of dynamical ensembles small by 
using a fast fermion discretization like the staggered one 
while having the advantages of a more computationally 
demanding regularization - like domain wall fermions or 
overlap fermions - in the valence sector. 

An alternative to the chiral expansion for describing 
the pion mass dependence of any hadronic observable is 
a Taylor expansion around a finite pion mass. In con- 
trast to the chiral expansion, the Taylor expansion is 
performed around a nonsingular point and has a finite 
radius of convergence. Typically this convergence radius 
is given by the distance of the expansion point Af^g to the 
chiral limit. Usually, the expansion is performed in pow- 



reweighted one ensemble to the physical point directly 



while the RBC-UKQCD collaboration (Aoki et al. 2011 ) 



followed a mixed strategy where the ensembles were first 
reweighted to the physical strange quark mass and a sub- 
sequent extrapolation to the physical pion mass was per- 
formed. 



The general idea behind reweighting (Ferrenberg and 



Swendsen 19881 is to reuse an ensemble produced with 
a certain set of parameters po = {Poi "^lol to obtain pre- 
dictions with a different set p = {/3, Wi}. As discussed in 
sect. |II.C[ gauge configurations U €U with the original 
set of parameters po in the action are produced according 
to the weight 



u;(po; C^) Of TT det M{miQ; U) x e 



-SgWo-.U) 



(149) 



ers of the pseudoscalar mass square Af^ 
in 



such that expectation values of observables ( 24 1 may be 



-"'^iTo resulting formed by just summing them over gauge configurations 



Mx=a + b{M^-Ml)- 



-cMl + d{Ml-Ml)^ (148) 



{o\ 



E 



ueu 



1 



(150) 



but the chiral behavior of baryon masses can also been 
fairly successfully described - at least within the statisti- 
cal accuracy of current data sets - by a linear expansion 



in M^ (Walker-Loud et al. 2009). Optimal convergence 



For a new set of parameters p, one may in principle 
circumvent the generation of a new ensemble with the 
weight 



is achieved in principle by placing the expansion point 
at the middle of the interval spanned by all simulation 



w{p; C/) ex TT det M{mi; U) x e 



■SG{fi;U) 



(151) 



points and the physical point ( Durr et al. 2008 Lellouch 



20091). Note however that from a practical perspective ^y reusing the old ensemble generated with the weight 



the choice of the expansion point M^ does not play a 
role in the fit itself as a redefinition of Af^ may be ab- 
sorbed by redefining the lower order fit coefHcients a and 



w{pq; U) from ( 149 1 and putting the ratio of weights into 
the observable 



b of ( 148 1 



(0)p- 



One may also try to fit ratios (Durr et al. 2008) or 



differences (Bietenholz et al. 2010b I of baryon masses in 






(152) 



Although ( 152 ) would in principle allow for a com- 



order to cancel common contributions and obtain a more 
regular chiral behavior. Further it is possible to study 
SU(3) breaking effects in baryon multiplets in the 1/Nc 
expansion (Jenkins et al. 2010), which offers an alterna- 



bined reweighting in both the coupling and the masses, a 
reweighting was carried out in the quark masses only in 



(Aoki et al. 2010 2011 1. For this purpose it is necessary 



five way of describing the chiral behavior of baryon mass 
multiplets. 

An alternative to extrapolating or interpolating results 
to the physical point is tuning the bare parameters of 
the lattice theory such that the physical point is directly 
reached, i.e. that the dimensionless combinations of con- 
tinuum variables mentioned above assume their physical 
value on the lattice. While recent advances in lattice dis- 
cretizations, algorithms and computer technology have 
made such an approach possible in principle, the compu- 
tational overhead that is associated with the parameter 
tuning is still large and the physical point is generally 
only reached within the precision of the tuning proce- 
dure. 



to compute ratios of fermion determinants at different 
quark masses. As it is prohibitively expensive to com- 
pute them exactly, stochastic methods were applied. 
Of course the reweighting method has its limitations. 



As one can see from ( 152 ), reweighting does exponentially 



enhance or suppress the weight of individual configura- 
tions in an ensemble with an exponent that is extrinsic, 
i.e. contains an explicit volume factor. As it is crucial 
for any observable to be computed on the relevant subset 
of configurations for the specific parameters used, these 
exponential factors should not be so large as to allow for 
one configuration to dominate the expectation value en- 
tirely. This in turn limits the allowed range in parameter 
space that one may reach safely with reweighting depend- 
ing on the original set of parameters po and the volume. 
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Within this safe range, the relative suppression of some 
configurations is that of effectively decreasing the statis- 
tics. As the relative weight factors are explicitly com- 



puted (Csikor et al. 2004), one has a very good handle 



on these effects. Note also that one could in principle by- 
pass some of the negative effects with a multihistogram 



technique (Ferrenberg and Swendsen 1989) 



Fig. [TT] provides an overview of currently used lat- 
tice QCD ensembles with respect to their position in the 
y'2M^ — M^ vs. Afjr plane. Note that in leading order 
of the chiral expansion the axes are proportional to \rm 
and y/fnTa. The location of the physical point also indi- 
cated. As fig.jTTJindicates, the physical point has already 
been reached. 



B. Continuum extrapolation 

The removal of the cutoff, also known as continuum 
extrapolation, is an unavoidable part of any lattice cal- 
culation that wants to make a statement about the under- 
lying fundamental continuum theory. The severity of the 
continuum extrapolation however depends very strongly 
on both the action used and the combination of scale 
setting observable and measured observable. 

As discussed in sect. |II.E[ the simplest gluonic action 
does already have cutoff effects of 0{a}^ that can be im- 
proved to at least 0{g^a'^) through various techniques. 
In contrast, the scaling behavior of the various fermion 
actions is typically not as good. Formally, staggered 
fermions and twisted mass fermions at maximal twist as 
well as exactly chiral fermions show 0{a^) continuum 
scaling while Wilson type fermions generically start with 
0{a) scaling. Only improved staggered fermions such as 
asqtad have a leading scaling behavior of 0{g^a?). There 
are several caveats to this statement however. 

For twisted mass fermions, 0(0?) scaling is realized 



at maximal renormalized twist (Aoki and Bar 2004 



Frezzotti and Rossi 2004a), which requires the tuning 
of one additional parameter. This tuning is routinely 
done as part of any twisted mass calculation (see e.g. 
(Baron et al. 2010b)). Since this tuning has a typical 
accuracy on the few percent level, it is expected that 
the 0{a?) terms are numerically dominant. In addition, 
twisted mass calculations often employ a doublet of va- 
lence fermions with oposite Wilson parameter to cancel 
remnant 0{a) effects. 

Similarly 0{a^) scaling is only strictly realized for chi- 
rally symmetric fermions if the chiral symmetry is exact. 
Fermion formulations that incorporate an inexact chiral 
symmetry, such as domain wall fermions at a finite fifth 
dimension, do formally have a remaining 0{g^"'a) scaling 
behavior. The smallness of the residual mass and other 
numerical evidence ( Aoki et al. 2011| ) however suggest 
that, similar to the twisted mass case, the 0{a^) term is 
dominant although it is formally subleading. 



Wilson type fermions on the other hand are typically 
Symanzik improved by the addition of a Sheikholeslami- 
Wohlert (clover) term (102). At tree level (csw = 1), this 
results in an 0{g^a) scaling of on-shell observables while 
with a suitable nonperturbative tuning one can in prin- 
ciple obtain 0{a'^) scaling. In addition, there is numer- 
ical evidence (Durr et al. 



2009 2011b Hoffmann et al. 



2007 Kurth et al. 2010b that the scaling behavior of 



clover fermions is substantially improved by gauge link 
smearing, which is commonly used today. 

Apart from the action used, the continuum scaling is 
also largely dependent on the observables considered. As 
discussed in sect. |IV.A] all observables in lattice QCD are 
dimensionless quantities and in order to extract dimen- 
sionful quantities such as baryon masses a scale setting 
observable is needed. The scaling is of course affected 
by the choice of scale setting variable. For baryons and 
vector mesons good scaling is observed when choosing a 
stable light baryon mass as the scale setting observable 



(Alexandrou ei aL 2009 Durr ei aL 2008) 



Some care has to be taken about the size of the scaling 
window. While generally scaling is not expected to set 
in for lattice spacings coarser than a ~ 0.1 — 0.15 fm, it 
has been observed ( Antonio et al.\ 2008 Bazavov et al. 



2010c 



et al. 



Del Debbio et a/.| |2002[ |Luscher[ [2010| |Schaefer| 



2011 ) that for fine lattices the autocorrelation time 



of the topological charge is rapidly growing. It therefore 
seems to be prohibitively expensive with current algo- 
rithms to obtain a sufficiently large and statistically in- 
dependent ensemble of configurations for lattice spacings 
finer than a ~ 0.05 fm. 

Generally, for the observables considered in this re- 
view continuum scaling is rather mild and not a leading 
source of systematic error. Fig. 12 gives an overview of 
the lattice spacing a vs. A/^ for currently used lattice 
ensembles. 



C. Finite Volume effects 

Besides reaching the physical point and removing the 
cutoff the third step that generically has to be taken in 
order to make physical predictions is the extrapolation to 
infinite volume. As is the case for the continuum limit, 
the infinite volume limit can never be reached and an ex- 
trapolation in the volume is in principle unavoidable. For 
most observables however the leading finite volume cor- 
rections are exponentially small in the box size and not 
polynomially and can therefore be made sufficiently small 

| l986a|). 



in practice by increasing the volume (Luscher 



These finite volume effects are discussed in sect. lIV.C.ll 
Resonant states on the other hand are embedded into 
a continuum of scattering states at infinite volume. In 
finite volume these levels become discrete and carry a 
strong volume dependence. Consequently the leading fi- 
nite volume effects on resonant states are of a different 
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I Baron erar||2010bl ) , MILC'IO (|Bazavov eTa l. 2010a'|, 



Data points are taken from the following references: ETMC'10(2+1+1) 



QCDSF-UKQCD'IO (iBietenholz et al.\ |2010a|), BMWc'08 purr et ~al 



2008D, B MWc'lO ( |Durr et aq|2011b|l, PACS-CS'0 9 



(Aoki et aZ.[ [2009a, 2010|, RBC-UKQCD'IO (|Aoki et a^|2011||Mawhinney[|20iol), JLQCD/TWQCD'09 ^sfoaki et fflZ.[p009|, 



HSC'IO ( |Lin et al.^ |2009| and all ensembles are from Nf = 2 + 1 simulations except explicitly noted otherwise. For staggered 
respectively twisted mass ensembles, the Goldstone respectively charged pion masses are plotted 



origin and are discussed separately in sect. |IV.C.2 

Finally we would like to mention that fixing the global 
topological charge in QCD is a restriction that becomes 
irrelevant in the infinite volume limit, too. For this rea- 
son lattice QCD calculations in a fixed topological sector 
may be viewed as introducing an additional third type of 



finite volume corrections ( Aoki et al. 2007 Brower et al. 



2003 ). Since at the time of this writing this technique has 



not been used in any work on light hadron spectroscopy 
we will not discuss it any further. 



1. Finite volume effects for stable particles 

In an interacting field theory, the properties of a par- 
ticle in a finite box are affected by mirror charge effects. 
For hadron spectroscopy this entails that all hadron 
masses in a finite box deviate from their infinite volume 
value with a leading contribution originating from the 
pion warping around one spatial lattice dimensionj^ A 



generic expectation for the finite volume correction to 
any hadron mass M in an L'^ x T box is thereforq^ 



Afoo 



(153) 



As (Luscher 1986a) demonstrated, there is a relation 



between the euclidean finite volume mass correction of 
a hadron P and the forward ttP scattering amplitude 
in Minkowski space. Concentrating on the case where a 
single propagator receives finite volume corrections, he 
obtained an explicit expression for the leading term in 
an expansion for asymptotically large L. Using an alter- 
native approach, (Gasser and Leutwyler 1987a|b 1988) 
incorporated finite volume effects into chiral perturbation 
theory. They demonstrated that the finite volume affects 
only the propagators and that it can be accounted for by 
simply replacing the momentum integration by a sum- 
mation over the allowed discrete momenta pi = 27:111.1/ L. 



Expanding the relation of (Luscher 1986a) to include 



subleading terms in asymptotic L and using xPT input 



^^ Alternatively in the momentum space view these effects may be 
considered as consequences of the discreteness of the momenta 
in a finite box. 



For the case of smaller volumes see also JFukugita et al.\ |1992| . 
They argue that the dominant (polynomial) finite size effect is 
due to the truncation of a hadrons wave function. 
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for the scattering amplitudes, (Colangelo and Durr 2004 2. Finite volume effects for unstable particles 



Colangelo et al. 2005 1 have combined the two approaches 
mentioned above for the case of pseudoscalar mesons. A 
similar expansion for baryons has also been pioneered 



(Colangelo et al. 2010) 



From a practical point of view these results imply that 
there is a safe asymptotic region of relatively large lattice 
volumes where these finite size effects are exponentially 
small and in addition can be systematically corrected for. 
As a rule of thumb for lattice computations with pion 
masses above ^ 300 MeV, lattices with M^L > 4 are 
considered safe while those with mT^L < 3 are widely 
affected by finite volume corrections. For a more quan- 
titative statement, fig. 13 shows a plot of box size L vs. 



pion mass Mjr where regions are identified that accord- 



ing to (Colangelo et al. 2005) imply the finite volume 



effect on the pion mass to be < 1%, < 0.3% and < 0.1% 
respectively. On top of these regions parameters of cur- 
rent or recent lattice computations are superimposed. As 
one can see current lattices are typically large enough 
to have percent level or smaller finite volume correc- 
tions on the pion mass. Note however, that corrections 



to baryon masses can be substantially larger ( Colangelo 



et al. 2010) 



Finite volume corrections are not always just expo- 
nentially small at large L as discussed in sect. |lV.C.l| In 
case where one is interested in extracting the mass of a 
resonant state that in infinite volume is embedded into 
a continuous spectrum of scattering states finite volume 
effects are more complicated. For illustration we start 
by considering the hypothetical case where there is no 
coupling between the resonance (which we will refer to 
as "heavy state" in this paragraph) and the scattering 
states. In a finite box of size L, the spectrum in the 
center of mass frame consists of two particle states with 
energy 



Af2 



fc2 



Af| 



fc2 



(154) 



where ki = 2ni7r/L and Mi, M2 are the finite volume 



masses of the lighter particles (cf. sect. IV.C.l) and, in 
addition, of the state of the heavy particle with finite 
volume mass Mx- As we increase L, the energy of any 
one of the two particle states decreases and eventually 
becomes smaller than the energy Mx of X. An analo- 
gous phenomenon can occur when we fix L but reduce 
the quark mass since the energy of the two light particles 
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changes more than Mx- In the presence of interactions, 
this level crossing disappears and, due to the mixing of 
the heavy state and the scattering state, an avoided level 
crossing phenomenon is observed. Such mass shifts due 
to avoided level crossing can distort the chiral extrapola- 
tion of hadron masses to the physical pion mass. 

The nterature (|Durr et al] |2008j |Luscher[ |1986b 
1991a|b| 



Rummukainen and Gottheb 1995) provides a 



conceptually satisfactory basis to study resonances in lat- 
tice QCD: each measured energy corresponds to a mo- 
mentum, |k|, which is a solution of a complicated non- 
linear equation. We follow (Luscher 1991a) where the 



yO-resonance was taken as an example and it was pointed 
out that other resonances can be treated in the same way 
without additional difficulties. The p-rcsonance decays 
almost exclusively into two pions. The absolute value of 
the pion momentum is denoted by k = |k|. The total 
energy of the scattered particles is 



W^2^/Mf+~k^ 



(155) 



in the center of mass frame. The tttt scattering phase 
Sii{k) in the isospin / = 1, spin J = I channel passes 
through 7r/2 at the resonance energy, which correspond 



to a pion momentum k equal to 



M2 



In the effective range formula 

(k^/W) -cottJii =a + bk^ 



this behavior implies 



'bki 






(156) 



(157) 



(158) 



where Tp is the decay width the resonance (which can be 
parametrized by an effective coupling between the pions 
and the p). The basic result of ( Luscher] 1991b) is that 
the finite- volume energy spectrum is still given by ( |155 ) 
but with k being a solution of a complicated non-linear 
equation, which involves the tttt scattering phase Sii{k) 
in the isospin / = 1, spin J — 1 channel and reads 



nvr- Sii{k) = 4){q) 



(159) 



Here k is in the range < fc < ^/ZM^r^ n is an integer, 
q = kL/(2Tr) and (j){q) is a known kinematical function 
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FIG. 14 Plot of the finite volume energy levels vs. box size 
1/ according to ( Luscher 1991b I. The full lines roughly corre- 



spond to the p—TTTv system, the dashed lines show the behavior 
in the case of a much smaller coupling gp. 



which can be evaluated numerically. In the limit of small 
q, (j){q) ex q^ and ^(q) w irq^ for q > 0.1 to a good ap- 
proximation. Solving the above equation leads to energy 
levels for different volumes and pion masses. 



Fig. [14] il- 
lustrates these solutions as a function of the box size. 

Thus, the spectrum is determined by the box length L, 
the infinite volume masses of the resonance Mx and the 
tvi^o decay products Mi and M2 and one parameter, gx, 
which describes the effective coupling of the resonance to 
the two decay products and is thus directly related to the 
width of the resonance. 

In infinite volume the resonance manifests itself as an 
increased state density in a continuous spectrum. Identi- 
fying the infinite volume resonance on a lattice at a given 
finite volume is therefore not straightforward - generically 
it is not possible to identify the resonance with a single 
energy level. Scanning over different system lengths it is 
however possible to identify energies with an increased 
probability of finding a statej^ as seen in fig. 14 This 
property may in principle be used to identify a resonance 



in a lattice calculation (Bernard et al. 2008b 2011 Giu 



dice et al. 2010). An alternative method based on fi- 



nite time correlators has also recently been suggested by 



(Meissner et al. 2011) 



Although conceptually clear, the treatment of resonant 
states in a region where they are not the ground state 
faces the huge challenge of reliably extracting the ground 
state as well as a number of excited states. One therefore 
often uses the assumption that an operator which does 



^'^ Note that depending on the specific volume this might be either 
the ground state or any of the excited states (cf. fig.|14|. 



not mirror the valence quark structure of a scattering 
state will almost exclusively couple to the resonance for 
extracting directly the desired resonance level. Recent 



studies (Engel et al. 2010 Lin et al. 2009) provide some 



evidence for the validity of this assumption. 



D. Electromagnetic effects and isospin breaking 

Including dynamical fermions into a lattice QCD cal- 
culation is numerically expensive due to the occurrence 



of the fermion determinant in ( 24 1 . For this reason, early 



lattice calculations have been performed ignoring the ef- 
fect of dynamical fermions all together (quenched approx- 
imation) . With improved algorithmic understanding and 
the increase in available computer power the inclusion of 
some dynamical fermion effects has become possible. The 
first step is usually the inclusion of a degenerate pair of 
light quarks followed by the inclusion of a single strange 
quark. 

Although the physical u and d quarks are far from de- 
generate mu/md ~ 0.56 (Nakamura et al. 2010), both 



masses are much smaller than the QCD scale Aqcd and 
therefore can be treated as a pair of degenerate light 
quarks with mass rh = {rriu + md)/'i to a very good ap- 
proximation. In addition, quarks are electrically charged 
and a full understanding of the experimentally observed 
hadron spectrum therefore necessarily includes QED ef- 
fects, too. For most observables related to hadron spec- 
troscopy, these are however subdominant as can be read- 
ily seen by the smallness of the coupling constant a-EM 
relative to the QCD coupling constant a. Effects of other 
interactions or of heavier quarks are negligible for light 
hadron spectroscopy within the currently attainable pre- 
cision. 

The relatively largest electromagnetic/isospin breaking 
effects can be observed in the pious and also for kaons the 
effect is still at the percent level. Since both pion and 
kaon masses are often used to define the physical point, 
it is necessary to define a properly isospin averaged pion 
and kaon mass as an input to lattice QCD calculations 
with a degenerate pair of light quarks. We denote the 
physical mass splittings in the pion/kaon sector as A^ = 



M^ 



ML and A^ = M 



/^± — Mj^o . Calling the pion 
and kaon mass in pure, isospin averaged QCD Mt^ and 
Mk respectively, we can write the pion and kaon masses 
in full QCD-I-QED as 



Mi 






and 



ll 



lk+ 



Tk+ = Mi 






(160) 



(161) 



where the Ix are the isospin and the F^; are the QED cor- 
rections to the squared mass of the particle x. In order to 
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obtain these, we start by noting that upon interchanging 
u and d in the pion sector 7r+ and tt~ are interchanged 
while for kaons a K^ goes over into a K'^ and vice versa. 
Disregarding QED effects, we therefore see that both 
M^4. and {M^^ + M^g) /2 may only contain even pow- 
ers when expanding in the isospin breaking parameter 
{rriu — rrid) and are therefore /^+ ex [niu — rrid)^ and 
Iko+Ik+ c< {rriu—md)'^ are free of leading (linear) isospin 
breaking effects. Moreover, as (Gasser and Leutwyler 
1984) have demonstrated, 1^^+ ex (m„ — miiY{m^ + rrid) 
leading to an even stronger suppression. Within cur- 
rently attainable precision, isospin breaking corrections 
to both Ml^ and {MJ.+ + Af^o) /2 are therefore negli- 
gible and we can assume I^+ = I^o + 1^+ = 0. 

Regarding the isospin correctio n to the charged pion, 
1^0 — 1^+ ex {rriu — md)^ itself and (Gasser and Leutwyler 



has been difficult to reach the physical point even in the 
isospin limit with two degenerate light quarks and since 
at least for the pion QED effects are substantially larger, 
it is these QED effects that have received most attention 
in the literature to date. 

Regularizing QED on the lattice poses a very differ- 
ent set of problems than regularizing QCD. There is a 
straightforward way of including the leading QED correc- 
tions into QCD calculations that was first employed by 



(Duncan et al. 1996) to obtain an estimate for the pion 



mass splitting. Since QED is an abelian gauge theory, 
the electromagnetic gauge field is trivial when ignoring 
the electrical charge of the sea quarks. In this partially 
quenched approximation one can therefore sample free 
QED C/(l) fields independently of the QCD SU{3) fields. 
This is most conveniently achieved by not sampling the 



1985) have found a parameter free expression at NLO in parallel transports Uf^{x) 



Ax)a directly but by in- 



terms of physical meson masses that yields 



/. 



IttO 



A. 



0.04 



(162) 



stead fixing the gauge and sampling the underlying gauge 
fields Afj_{x) in which the Maxwell equations are linear. 
The corresponding action 



indicating that the bulk of the pion mass splitting is due 
to electromagnetic effects. 

Regarding electromagnetic effects, according to 
(Dashenj 1969) the picture in the leading order of the 
SU{3) chiral expansion is that while neutral pseudo- 
Goldstone masses stay unaffected by electromagnetic cor- 
rections, i.e. Tt^o = F/fo = 0, the difference of the 
square of the charged masses receive the same correction 
r^+ — Tx+ ■ The absence of electromagnetic corrections 
in the 7r° and to a lesser extent the iC" mass is further jus- 



'QED 






(d^A.ix) ~ d,A^{x)Y 



(164) 



tified by (Das et al. 1967) who demonstrated that these 
corrections vanish in the massless limit m„ = m^ ~ 
respectively rriu = rrid = trig = 0. On the other hand, a 



range of model calculations (Bijnens and Prades 1997 



Donoghue and Perez 1997) that are partly based on the 



inclusion of QED effects into chiral perturbation the- 



ory ( ,Urech 1995) and dispersive calculations based on 
the 77 — >■ Stt decay (Anisovich and Leutwyler 1996| Bij- 



with the forward difference operator d^ decouples in 
Fourier space and gauge field configuration with the 
proper weight e~ '^^'° may thus be produced by simply 
producing each Fourier component with a proper random 
weight. The only further restriction is the vanishing of 
the p — component due to the compact space provided 
by the finite lattice. 

In a calculation with dynamical sea quarks, the elec- 
tromagnetic corrections to the light quark masses that 
is introduced in the valence sector leads to a mismatch 
of sea and valence quark masses. In order to minimize 
these unitarity violating effects, one can retune the va- 
lence light quark masses such that after the inclusion of 
quenched QED effects both rriu and rud have the same 



nens and Ghorbani 2007[ Colangelo et al., 2009 ; Ditsche 



et al. 2009 Kambor et al. 1996 Leutwyler, ,1996) sug 



value as they had in pure QCD (PorteUi et al. 2010). 



In the quenched approximation calculation of (Dun 



gest that there are noticeable corrections to the other 
parts of Dashen's theorem. A recent world average of 



can et al. 1996) a large violation of Dashen's theorem 
was found corresponding to e ^ 0.5. In later work 
with Nf = 2 dynamical flavors of domain wall fermions. 



(Blum et al. 



2007) have reported a somewhat larger 
these corrections was given in ( [Colangelo et al.\ |2011[ ) as value whilem the most recent update Nf =2 



(r 



RO 



K+ 



et al. 



1 (Blum 



2010 1 foimd a value compatible with the origi- 
0.7(5) (163) nal estimate of (Duncan et al. 1996). Preliminary re- 



Assuming Dashen's theorem to hold we would thus get 
for the QED corrected, isospin averaged pion and kaon 
masses M^ ~ 134.8 MeV and Mk c^ 495 MeV. Includ- 



sults are also available from the MILC collaboration us- 
ing Nf = 2 + 1 flavors of staggered fermions (Basak et al. 



ing the correction ( 163 ) while keeping the reasonable as- 2010 ) 



2008) and the Budapest-Marseille- Wuppertal collabora- 
tion with Nf = 2-1-1 Wilson fermions (Portelli et al. 



sumption Tj^o = T^o — 0, the kaon mass is shifted by 
less than 1% to Mk — 494.6 MeV while the pion mass is 
unaffected. 

In order to go beyond these results, we need to treat 
QED and isospin breaking effects on the lattice. Since it 



The general picture that emerges is that the correc- 
tions to Dashen's theorem that are parameterized in e are 
in agreement with the phenomenological determinations. 



Taking lattice determinations into account, (Colangelo 



et al. 2011) have combined recent results on QED and 



37 



isospin splitting effects into a world average of the QED 
corrected, isospin averaged pion and kaon masses. For 
the individual symmetry breaking parameters they find 



e = 0.7(5) e,„ = 0.04(2) 

r^o = 0.07(7)A^ Tko = 0.3(3)A^ 



(165) 



and consequently M^^ = 134.8(3) MeV and Mk = 
494.2(5) MeV, which agree within error with the values 
quoted above that were obtained under the assumption 
that Dashen's theorem holds. 

Isospin splitting effects in the baryon spectrum are 
less dramatic than for the pseudoscalar mesons. As an 
example, the mass splitting in the nucleon system is 



Af„ 



1.2933321(4) MeV (Nakamura et al. 2010) 



and an effective theory estimate of the electromagnetic 
contribution turns out to be negative {Mn — Mp)^^^^ = 



-0.76(30) MeV (Gasser and Leutwyler 1982 1. If one 



is interested in isospin averaged baryon masses only, a 
straight isospin average is therefore sufficient for the ac- 
curacies that are presently obtainable in lattice calcula- 
tions. 

A first dedicated lattice study of the nucleon mass 



splitting was carried out by (Beane et al. 2007). They 



used a hybrid setup with domain wall valence on Nf = 
2-1-1 staggered sea quarks with pion masses in the range of 
~ 290 - 350 MeV. A single lattice spacing a ~ 0.125 fm 
was used and an extrapolation to the physical point was 
carried out in the framework of NLO partially quenched 
heavy baryon xPT. Using the experimental value of 
Ma — Mjv as an input they obtain for the isospin part of 
the mass splitting (M„ - Mp)^^^ = 2.26(57) (43) MeV. 
In addition to the isospin part of the nucleon mass 



diference, (Blum et al. 2010) have also calculated the 
QED part. They use Nj = 2 + \ partially quenched 
domain wall fermions with pion masses in the range 
of ~ 250 — 400 MeV at a single lattice spacing a ^ 
0.11 fm. Two volumes were used to estimate finite 
size effects in the final result and an extrapolation to 
the physical point was performed using NLO partially 
quenched heavy baryon xPT. They quote the final results 
(M„-Mp)qcd = 2.24(12) MeV, (M„ - Mp)^^^ - 
-0.383(68) MeV and M„ - Mp = 1.86(14)(47) MeV 
where the first error is statistical and the second part 
of the systematic error. 



V. LATTICE RESULTS 

In this section we discuss lattice results on the light 
hadron spectrum. Historically, the first results were 
from the quenched approximation that is discussed in 
sect. |V.A[ The inclusion of dynamical fermions was pi- 
oneered with heavy, degenerate quarks (sect. |V.B | be- 
fore it developed into the study of theories with non- 
degenerate light and strange quarks that we review in 



sect. V.C While sect. V.A and |V.B| are now of mainly his- 
torical interest, the three-fiavor (and coming four- flavor) 
dynamical calculations are the definitive modern calcu- 
lations. 

Our review of lattice results does not include glueballs. 



For recent reviews on this topic see e.g. (Klempt and 



Zaitsev[ [20071 |Mathieu eTal] [20091 |McNeile|l2009| |Teper 
1998[). 



A. Results in the Quenched approximation 

Although the quantitative understanding of the light 
hadron spectrum is an obvious and essential check for any 
candidate theory of the strong interaction, it took several 



years from the original proposal of (Wilson 1974) that 



contained the lattice discretization of gauge theories and 
the strong coupling picture of quark confinement to the 



first numerical studies of the hadron spectrum (Fucito 



1981 1983^ Hasenfratz et al. 1982 Marinari et al.\ 1981a 



et al.\ [19821 [Fukugita et al.\ [19841 [Hamber and Parisi 



Martinelh et al. 1982 Weingarten 1983). Due to the 



lack of viable dynamical fermion algorithms and com- 
puter power, these pioneering studies were carried out in 
the quenched approximation sometimes with an SU(2) 
gauge group and even further discrete truncations. Lat- 
tices had a typical size of 6^ x 12 and 0(10) gauge con- 
figurations were generated with the Wilson gauge action. 
Naive or plain Wilson fermion actions were typically used 
to extract hadron masses and physical point predictions 
were obtained by linear extrapolation of either squares 
of meson masses or baryon masses. A first world average 



of these pioneering results was given by (Creutz et al. 



1983) 



rup = 800(100) MeV 
TOao = 950(150) MeV 
TOai = 1100(150) MeV 

THp = 1000(150) MeV 
ruA = 1300(150) MeV 



(166) 



From a modern perspective, these results should be 
viewed with some caution as these calculations were 
clearly exploratory and pioneering. The computer power 
of the times was not sufficient to properly clarify many 
systematic effects. As an example, the inverse lattice 
spacing of SU{3) gauge theory with the Wilson gauge 
action at /3 — 6.0 used by (Marinari et al. 1981a) was 



a — 1.12 GeV whereas modern determinations from 
various observables agree that it is a~^ ~ 2.1 — 2.3 GeV 



(Aoki et al. 



2000 



Durr et al. 2007 Giusti et al. 



mer 2002 Otto and Stack 1984 1. 



It was quickly realized ( [Bernard et al 
T9831 [Gupta and Patell [1983 



et al. 



Montvay 



2001 



Gutbrod et al.[ T983[ Lipps et al. 1983| Necco and Som 



1983 




Politzer 



1983 Bowler 



asenfratz and 



1984) that physical volumes 
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were not big enough and that one should use larger time 
extents in order to safely extract a ground state signal. 
In the following years, quenched calculations with unim- 
proved Wilson and staggered fermions on Wilson gauge 
action were pushed to larger lattice volumes and higher 



statistics (Billoire et aL. 1984 1985 Bowler et al. 1984 



19851 [Gilchrist et aL[ |1984| |Itoh eraLH1986a|b[ JKonig 



et ai.[ |1984[ [Kunszt and Montvay| |1984[ [Lipps et al. 
1983) where lattices were often doubled in time direc- 
tion in order to obtain a clean signal. With gauge cou- 
plings typically /3 ^ 5.7 — 6 and spatial lattice extents 
typically 10—16 lattice units and time extents typically 
twice as much, a qualitatively consistent picture of the 
hadron masses started to emerge although large system- 
atic effects were present that could not clearly be iden- 
tified yet. (For reviews of this generation of results see 
( Hasenfratz and Hasenfratzj 1985| Montvay 1987)). In 
particular, the ratio of the nucleon mass to the p mass, 
which experimentally is M^/Mp ~ 1.21 turned out to 
be consistently too high M-^ jMp > 1.6, which is even 
larger than the static quark limit Mn/Mp = 1.5. An- 
other stumbling block for these early calculations was 
the absence of sufficient splitting between the masses of 
the nucleon and the A. 

During the following years, the focus shifted slightly to- 
wards inclusion of sea quark effects with steady progress 



in quenched spectroscopy (Bacilieri et al. 1988 1989 



Fukugita et al. 1988 Gupta et al. 1987) until the ffist 



precision calculations of the quenched light hadron spec 
trum emerged in the early 1990's ([AUton et al. 1992 



T994^ 'Bacilieri et al.\ 1990; Bitar et al., 1992; Butler et al. 
1993, ,1994, iCabasino et aL,|1991a,b,|Daniel et aLH1992 



clair 



1993). 



Guagnelli et al. 1992 Gupta et al. 1991 Kim and Sin- 



Among these, the first landmark precision calculation 
of the quenched light hadron spectrum was carried out 



by the GFll collaboration (Butler et al. 1993 1994 1. 



Wilson fermions were used on a Wilson gauge action at 
three different values of the lattice spacing in the range 
a ^ 0.07 — 0.14 fm. Propagators were extracted using 
point and gaussian smeared sources at different quark 
masses corresponding to M-j^/Mp > 0.5 i.e. with pion 
masses M^. > 400 MeV. Lattice volumes in the range 
16"^ X 32 to 30^ X 32 X 40 were used corresponding to 
a spatial lattice extent of ^ 2.3fm at all three lattice 
spacings. At the coarsest lattice spacing, dedicated runs 
at larger and smaller volumes were performed in order to 
extract the finite volume dependence of the result. They 
were used in the end to correct the physical predictions 
to infinite volumerj Considering degenerate quarks only, 
a linear relation was established between the degenerate 



^* See ( |Gottlieb[|l997| for a detailed discussion of the finite volume 
effects. 



ratio 


finite volume 


infinite volume 


observed 


ijiK'/mp 


1.149±0.010 


1.167±0.016 


i.ie4 


m^/mp 


1.297±(l.019 


1.333±().032 


1.327 


m^i/mp 


1.285±0.070 


1.219±().105 


1.222 


Am/mp 


1.867±0.046 


1.930±0.073 


2.047 


mf^lnip 


1.628±0.075 


1.5fJ5±0.111 


1.6(14 


m^^./nip 


1.813±0.051 


1.821±0.075 


1.803 


in^-jinp 


2.013±0.052 


2.063±0.067 


1.996 


ma/mp 


2.206±0.05S 


2.29S±0.098 


2.177 


A— /"V 


O.3()5±().0O8 


0.319±().012 


0.305 ± 0.018 
0.32(1 ± 0.007 



TABLE I Quenched lattice QCD prediction of the light 
hadron spectrum according to (jButler et al. 1993 19941. The 



table gives the ratio of various hadron masses and the QCD 
scale Ai-p4 to the mass of the p that was used to set the scale. 

MS '^ 

The label Am refers to the combination Am = ms+m-E—m,N. 



Observed values are experimental results from (Hikasa et al. 



i((!T 



1992|) except for the case of the QCD scale Ai-i- where they 



refer to two previous results from the literature ([Bali and 



SchillingI |1993a|b[ |E1-Khadra et aL\ |1992[ ). Note that some 



experimental values, notably the mass of the p, have been 
updated since (Nakamura et al. 20101. 



quark mass iriq and M^ while for all other hadrons a fit of 
the form M = a+bniq described the data. Assuming that 
these linear relations extend to the non-degenerate case 
with two quarks of mass mi and m2, i.e. M^ oc {mi + 
TO2)^ and M = a + biirii + h2m2 the physical point was 
found using M^r, Mx and Mp input with the later used 
as scale setting observable. A continuum extrapolation 
linear in a was performed that turned out to be rather 
mild. 
Tablell 



ler et al. 



shows the resulting spectrum obtained by ( But 



1993 1994). Despite the many approximations 



used, the overall agreement with experiment is rather re- 
markable and at the < 10% level. 

Similarly sophisticated quenched analyses were soon 
after performed for the 77 — 77' system ( Kuramashi et al. 



1994) and for excited state mesons (Lacock and Michael 



1995 



Lacock et al. 1996). These calculations and de- 



tailed studies of systematic effects such as finite size 



(Aoki et al. 1994), excited state contaminations (jlwasaki 



et al. 1996) or quenched chiral logarithms and S'C/(3) 



splittings ( Bhattacharya et al. 1996[ |Kim and Sinclair] 
1995 ) revealed potential inconsistencies of the quenched 
approximation of up to 20%. On the other hand, results 



with improved Wilson actions (Collins et al. 1997 Ed- 



wards et al. 1998 Gockeler et al. 1997 1998b) and for 



staggered fermions that reached finer lattice spacings and 



smaller quark masses (Bernard et al. 1998a Kim and 
Ohta 2000[ ) indicated quenching effects that were less 
dramatic at 0(5%). In the case of the latter two results 
it was especially noted that a simple linear extrapolation 
in the light quark mass was no more sufficient. Several 
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FIG. 15 



Quenched lattice QCD prediction of tlie ligiit hadron 

For both 



2000 2003b I 



spectrum according to (Aoki et al 
data sets plotted Mp and M-^ were used to set the scale and 
the light quark mass. In order to set the strange mass, ei- 
ther Mk (filled circles) or M^ (open circles) are used. Plot 
reproduced with friendly permission of the CP-PACS collab- 
oration. 



xPT motivated fit forms were found to describe the chiral 
behavior of the p and nucleon masses but the coefficients 
were not found to be in agreement with quenched xPT 
expectations at aU. 

The accuracy of the quenched approximation was ad- 
dressed in the large scale calculation by the CP-PACS 

They used lat- 



2000 2003b) 



collaboration (Aoki et al 
tices of ^ 3 fm spatial extent at four values of the lat- 
tice spacing in the range a ~ 0.05 — 0.1 fm with quark 
masses corresponding to M^/A/p ^ 0.4 — 0.75. Both 
fermion and gauge action used were plain Wilson and 
non-degenerate quark masses were used to investigate 
splittings in the SU{3) multiplets. While pseudoscalar 
meson masses were found to have a chiral behavior com- 
patible with the quenched xPT expectations, the dis- 
crepancy in the vector meson and baryon sector found in 



the staggered results of (Bernard et al. 1998a Kim and 



Ohta 2000) was confirmed. For these masses, xPT moti- 



vated fits were used to extrapolate to the physical point. 
The physical point in the light quark mass was defined 
using Mtt and Mp and either Mk or M^ were used to 
define the physical strange massj^ The final result that 
has a precision of ^ 1 — 3% is displayed in fig. 15 



A statistically significant deviation from the experi- 
mentally observed spectrum was noted with discrepancies 



^^ Note that although the </> is a singlet meson, its disconnected 
part is usually disregarded in lattice studies. 



up to ~ 10%. This discrepancies however are particularly 
pronounced due to the choice of Mp as a scale setting ob- 
servable. Since the p does not decay in the quenched ap- 
proximation and therefore represents the ground state in 
the vector channel, it is in principle a viable scale setting 
variable from a pure lattice perspective. Nonetheless, the 
identification of the stable quenched ground state energy 
with the mass of a resonance with ^ 150 MeV experi- 
mental width is not optimal. On top of that, an accu- 
rate determination of p meson properties is a challenging 
experimental task. This is highlighted by the fact that 
the experimental value of Mp itself has moved by ^ 1% 
or more than ten standard deviations over the last two 



decades (Hikasa et al 1992 Nakamura et al. 2010) 



As ( Garden et al. 2000 ) noted, one can derive from the 



CP-PACS results predictions for hadron masses with the 
scale set by the nucleon mass instead of Mp. In this case, 
the maximum deviation from experiment turns out to be 
significantly lower at ~ 4% indicating that indeed the 
quenched approximation is substantially worse for reso- 
nance masses than for masses of hadrons that are stable 
within QCD. A confirmation of these results with some- 



what lower statistical accuracy was reported by (Bowler 



et al 2000) 



The CP-PACS calculation was one of the last large 
scale calculation aimed at a precision determination of 
the quenched ground state light hadron spectrum. As 
quenching effects had become statistically significant, the 
focus of efforts to get a quantitative confirmation of QCD 
reproducing the experimentally observed ground state 
light hadron spectrum moved towards inclusion of dy- 
namical fermion effects. Nonetheless, due to its relatively 
low numerical cost, the quenched approximation contin- 
ues to be used to this day as a testbed for new numeri- 
cal approaches and as a first step in studying computa- 
tionally demanding observables. In the following years, 
the quenched ground state hadron spectrum was used to 
check various variants of chirally symmetric fermion ac- 



tions ([Babich et al 2006 Galletly et a/. 2007 Gattringer 



et al 2004 1 or to develop improved ( Melnitchouk et al. 



2003) or anisotropic (Nemoto et al 2003) actions that 



were intended for studying excited hadrons. 

Following the CP-PACS determination of the ground 
state quenched hadron spectrum the attention in 
quenched hadron spectroscopy turned towards resonant 



and singlet states (Basak et al. 2007 Brommel et al. 



20041 IBurch et a/.[ |2006a|b| |Engel et al\ |2010| |Flem 



ing e^ aLj |2009| IGattringer et aZ.[ |2004[ [Gockeler et al 



2002[ [Lasscock et al. 



Mathur et aLl p05 



Melnitchouk et al 



et al. 



2005 



2002 



[2007; Lee and Weingarten 



2007j .McNeile and Michael 



2003 



Nemoto et al 2003 



1999 



2001 



Sasaki 



Wada et al. 2007). In particular, many 



groups reported on the splitting between the nucleon and 
its lightest negative parity partner, the A^*(1535). In all 
cases, a clear signal of the mass splitting between the 
nucleon ground state and the A^*(1535) could be seen 
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on the lattice. The sphtting is increasing as one lowers 
the light quark masses towards the physical point and 
in all cases is roughly consistent with the experimentally 
observed mass splitting. It is an interesting peculiar- 
ity however, that the lightest experimentally observed 
nucleon excited state is not the nucleon parity partner 
iV*(1535) but in fact the iV*(1440), the so-called Roper 
resonance which carries positive parity and J = 1/2 as 
the nucleon. With the exception of ([Mathur et al. |2005[ 



Sasaki et al. 2005 1 who employed Bayesian techniques to 
extract excited state information from a single channep^ 
however, the first positive parity excitation of the nucleon 
turned out to lie above the first negative nucleon state 
in all lattice calculations. A possible solution to this dis- 



crepancy has recently been proposed by ( Mahbub et al. 



2010a 2009b I who demonstrated that a mix of excited 



states enters typical interpolating operators. By using a 
large operator basis they could explicitly disentangle up 
to eight states and demonstrate a level crossing between 
the negative parity ground state and the first positive 
parity excitation for light quark masses that is consistent 



with the finding of (Mathur et al.\ 2005 Sasaki et al. 



2005 ) and the experimentally observed level ordering be- 



tween the iV*(1535) and iV*(1440) states. 

Another interesting case in excited state baryon spec- 
troscopy is the A(1405). In the quark model picture, 
it is the lightest negative parity partner of the A with 
a valence quark structure uds. It is however the light- 
est negative octet baryon - more than 100 MeV lighter 
than the lightest negative parity nucleon, the A^*(1520) 
- even though it does contain a strange quark. This is 
the most striking of many peculiar features that have 
given rise to a number of suggestions for a nontrivial 
structure of the A(1405) such as that of a NK hadronic 
molecule or a pentaquark state (for a recent review see 



(Klempt and Richard[ 2010)). (Melnitchouk et al. 2003 



Nemoto et al. 2003) have studied the negative parity 



A states in the quenched approximation with standard 
interpolating operators and found it impossible to repro- 
duce the A(1405) which they interpret as an indication 
for a nontrivial structure of the A(1405) that might not 
be properly reflected in the quenched approximation. In 
contrast, ( Burch e^ oH |2006a| ) found the A(1405) to be 
consistent with the negative parity octet state. 

With the exception of the A^*(1440) and the A(1405) 
that were discussed above, no qualitative tension between 
experiment and the above mentioned quenched excited 
baryon studies was found. For a more detailed review 



see (Leinweber et al. 2005). 



B. Results with degenerate dynamical quarks 



The spectrum calculation of the CP-PACS collabo- 



ration (Aoki et al. 2000 2003b I reached a numerical 



precision such that quenching effects could clearly be 
seen. In order to obtain a quantitative understanding 
of the ground state light hadron spectrum on the few 
percent level it is therefore necessary to include dynam- 
ical fermion effects into the lattice calculation. From 
the six fermion flavors in nature, the charm, bottom 
and top each have masses much larger than the QCD 
scale Aqcd- Their contribution to light hadron masses 
through quark loop effects is therefore believed to be 
negligible. Among the remaining three flavors, rriu/d <?C 
Aqcd while m^ ^ Aqcd- Consequently and because an 
even number of quark flavors is usually easier to imple- 



ment (cf. sect. II.C), the first attempts at unquenching 



lattice QCD calculations were performed with two de- 
generate quark flavors. In the case of staggered fermions 
calculations with four dynamical quark flavors are even 
easier due to the remnant doubling (cf. sect. [II.D.l I. 
In this section we review results obtained with (a usu- 
ally even number of) degenerateFH dynamical fermions. 
While this still represents an approximation that is ne- 
cessitated by the lack of proper computational resources, 
it is still a very significant step forward from the quenched 
approximation. 

Pioneering work on lattice hadron spectroscopy with 



dynamical fermions was done by (Langguth and Mont 



vay 



1984 ) where dynamical fermions were implemented 



using a strong coupling expansion of the determinant ra- 
tio. During the following years unquenching via the in- 
clusion of Nf = 2 and 4 dynamical Wilson and Nj = 2, 3 
and 4 staggered fermions was investigated by several 



groups (Billoire and Marinari 1987 Campostrini et al. 



T987| IDetar and Kogut| |1987a|b[ (Fucito et d. 



1986 



Fukugita et aZ l| 1987a';i986;'1987b [|Gottlieb et aZ.[|1987a 



Grady et al. 1988; Hamber^ ,1987 1. These early works 
demonstrated that the main effect of including dynam- 
ical fermions was a change in dependence of the lattice 
spacing on bare coupling constant j3. Apart form this 
effect, no clear sign of unquenching could be observed. 
In particular, the Mp/M^ ratio tended to stay constant 
or even increase. In studies with staggered fermions, the 
taste breaking effects were also observed to be rather se- 
vere. For a comprehensive review of these early studies 



see (Fukugita 19881. 



During the following years it became clear that with 
staggered fermions one could go to substantially lighter 
quark masses than with Wilson fermions ( Altmeyer et al\ 



See fSasaki a nd Sasaki] |2005| for a discussion of possibly large 
finite size effects for this technique 



■^^ We use the term non-degenerate here in the sense that no explicit 
term was added that does break the flavor symmetry. In the case 
of staggered quarks the unavoidable taste splitting is of course 
present. 
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19931 IBitar et al.\ '1990a' b[ [T992{ [Brown et aL\ [1991 



Fukugita etal\ [1993, .Gupta et al.\ [T991J [Patel et al. 



1989) in the sense that the M^r/Mp ratio attainable with 



Wilson fermions was hmited to M^/Mp > 0.7 while if 
taking the Hghtest pion one could go down to about 
half this number in the staggered case. Since reduc- 
ing the mass of the valence quarks only was substan- 
tially easier, some of these studies started exploring par- 
tially quenched setups, where the valence quark masses 
are varied independently of the sea quark masses, and 
even hybrid calculations with valence Wilson quarks on 
a dynamical staggered sea. None of these calculations 
however gave a clear signal for a Mp/M^ ratio that was 
substantially better than the ones obtained in contempo- 
rary quenched calculations. Although there was steady 



progress over the following few years ( Allton et al. 1999 



Bitar et al. 1994 Bicker et al. 1998), the focus of large 



scale calculations shifted more towards precision compu- 
tations in the quenched approximation. This was in large 
part due to the tremendous computational effort that 
was needed for dynamical fermion computations which 
exceeded the computer capabilities of that time. A first 
unquenched study of the rj — rj' mixing was performed 



by (McNeile and Michael 2000) which found a mixing 



angle oi 6 ^ —10° albeit without continuum and chiral 
extrapolation[^ 

These efforts culminated in the first large scale project 
to compute the light hadron spectrum in Nf = 2 QCD by 
the CP-PACS collaboration ([Ah Khan et al.\\2002L They 



used two degenerate flavors of mean field improved clover 
fermions on an Iwasaki gauge action. The strange valence 
quark was included in a quenched setup. Three relatively 
coarse lattice spacings in the range a ^ 0.11 — 0.22 fm 
were used with an approximately constant physical vol- 
ume L ~ 2.5 fm and T — 2L. Four sea and valence quark 
masses in a range corresponding to Mp/M^q ^ 0.6 — 0.8 
and an additional valence quark mass at AIp/Mj^ ^ 
0.5 were investigated. Point sources and exponentially 
smeared quark sources on a gauge fixed background were 
chosen for optimal plateau onset. Chiral extrapolation 
was performed by a combined fit to all partially quenched 
masses for each channel on a given sea quark mass. Vec- 
tor meson and baryon masses were extrapolated to the 
physical point using quadratic functions in the valence 
and sea M^ with certain restrictions on the quadratic 
terms. In the case of vector mesons, xPT motivated 
nonanalytic M| type terms were also used instead of M^ 
type terms to estimate the systematic error. Following 
the example of the quenched CP-PACS calculation dis- 
cussed in sect. V.A the physical light quark masses and 



the scale are defined via M^r and Mp while two options, 
Mk or M^, were used to set the strange quark mass. The 



continuum limit is obtained by linear extrapolation in a. 

The resulting light hadron spectrum is plotted in 
fig. [TBI Clearly the heavier baryon states are in good 
agreement with experiment while the lighter ones, espe- 
cially the nucleon and the A, seem to be systematically 
too high. This does not come as a big surprise though 
since the extrapolation to the physical point is substan- 
tially more severe for the baryons containing more light 
valence quarks. 

Similar efforts to that of the CP-PACS collaboration 
were reported by the UKQCD and JLQCD collaboration 
in ([Allton et a/.[|2002[[Aoki et a/.[[2003a[). The UKQCD 



collaboration worked at a single lattice spacing a ~ O.lfm 
that was set with tq. The range of sea quark masses was 
chosen such that M^/Afp ^ 0.55 — 0.9 and the spatial 
lattice extent was L ~ 1.7 fm. The JLQCD collabora- 
tion worked at one single lattice spacing a ~ 0.09 fm at 
a spatial extent L ^ 1.8 fm using clover fermions on a 
Wilson gauge action and the same range of sea masses 
KU/Mp - 0.6 - 0.8 than CP-PACS. The findings of both 
collaborations on the light hadron spectrum are in good 
agreement with the continuum CP-PACS results. 

The conclusion from these two large scale projects re- 
garding the feasibility of computations with light dynam- 
ical quarks was summarized in a plot that became to be 
known as the "Berlin waU plot" by ([Ukawal |2002[). He 



conjectured that the cost of dynamical fermion simula- 
tions would rise towards the chiral limit essentially as 
ex M^ effectively rendering any dynamical calculations 
with Wilson type fermions near the physical point im- 
possible in any foreseeable future without substantial al- 
gorithmic improvements. Due to their lower computa- 
tional cost staggered fermions were able to push further 
towards the chiral limit and did already do so with also a 



dynamical strange quark (see sect. V.C), but extracting 



especially baryonic states is less straightforward in this 



formulation (cf. sect. III.Al. For a more comprehensive 
review of these results see ( McNeile 2003 ) . 



Because the inclusion of a non-degenerate sea quark 
incurs little extra expense, from this point on, the main 
development in light hadron spectroscopy continued with 
the inclusion of a strange quark in addition to two de- 
generate light quarks. This is usually referred to as 
Nf = 2-1-1 and is discussed in sect. V.C The question of 



how much two-flavor calculations differ from experiment 
has not been answered as definitively as it has been for 
the quenched calculations discussed earlier. A number 
of Nf = 2 calculations of light hadron masses were still 
being performed however for a number of reasons such 
as algorithmic tests (Del Debbio et al.. 2007a|b ) or as 
a first step for formulations that allow even number of 
quark flavors only ( Alexandrou et al. 2009 ) . Similarly, 



See also i Lesk et al. 2003 1 



flavor singlet spectroscopy that typically requires sub- 
stantially more statistical precision has still been inves- 
tigated in Nf = 2 QCD ([Allton e^ "alj [2004[ [Hart et al. 



2006 Kunihiro et al. 2004 McNeile and Michael 2006 
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FIG. 16 Dynamical Nf = 2 lattice QCD prediction of the light baryon spectrum according to (Ali Khan et al. 2002'). The two 



figures display the continuum extrapolation of the ground state light octet and decuplet baryon masses in the case where the 
strange mass was set via Mk ■ The continuum extrapolation was performed using three lattice spacings that are displayed with 
solid circles. Open circles represent a fourth, finer lattice which was not included in the analysis due to a too small volume. 
Plot reproduced with friendly permission of the CP-PACS collaboration. 



McNeile et aL, 2001 2009 Prelovsek et al. 2004). For 



a recent comprehensive review of these results see (Mc 



Neile 20071. 



The ETM collaboration (Alexandrou et al. 2009) has 



published results for the ground state light baryon spec- 
trum with Nf = 2 twisted mass fermions ^ They used 
Nf — 2 twisted mass fermions at maximal renormalized 
twist on a tree level Symanzik improved gauge action. 
Two lattice spacings (a ^ 0.07 fm and a ~ 0.09 fm) 
were used with charged pion masses in the range 270 to 
500 MeVJ^ The lattice spacing was set via the nucleon 
mass and chiral extrapolations were performed with a 
variety of different ansatze. The valence strange quark 
mass is set by tuning the kaon mass to its physical 
value. The final result employs two different heavy bar- 
ion xPT ansatze (O(p^) respectively NLO SU{2)) for 
the extrapolation of baryons without respectively with 
valence strange quarks to the physical mass point. The 
continuum extrapolation was performed using a constant 
which was demonstrated to be sufficient at the given level 
of accuracy. Exponential finite volume corrections were 
taken into account in the final fit form. Resonant state 



^^ Ses also J Alexandrou et al.\ |2008[ l for some results with more 

lattice spacings and different scale setting. 
^* The isospin splitting of the pio ns is M.^ — M^ 

220 MeV)2 (iBaron et aZ. I |2010cl 



(150 



finite volume corrections were not performed but are be- 
lieved to be irrelevant in the region of parameter space 
covered by the simulations. Effects of the twisted mass 
isospin breaking were observed to be negligible except in 
the case of the S where they amounted to a 6% correc- 
tion. Their final result is displayed in fig. [17] shows good 
agreement with experiment at the level of precision of the 
calculation which is ~ 5%. The ETM collaboration also 



investigated the p-uj mass splitting and mixing ( McNeile 



et al. 2009) excluding electromagnetic effects. While a 



clear signal and qualitatively correct behavior was found, 
the quantitative understanding of the experimentally ob- 
served splitting remains a challenging task. 

Turning to excited states, the BGR collaboration has 
computed ground and excited state hadron spectra us- 
ing Nf — 2 single step stout smeared chirally improved 
fermions on a tadpole improved Liischer-Weisz gauge ac- 
tion at a singe lattice spacing a ~ 0.15 fm (Engel et al. 



2010 Prelovsek et al. 2010). Three pion masses in the 



range 320—530 MeV were used and the scale was set with 
tq . Gaussian smeared quark sources were used in combi- 
nation with a variational method based on three interpo- 
lating operators to extract the energy levels. A chiral ex- 
trapolation linear in M^ was performed and the strange 
quark was introduced in a partially quenched setup. The 
results for positive and negative baryon states are plot- 
ted in fig. [18] A good signal for the ground state was 
found but excited and scattering state signals were gen- 
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by the BGR collaboration with Nf = 2 chirally improved 
fermions. The plot is reproduced from (Engel et al. 20101 
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erally weak. Some evidence was also presented that the 
a and k resonances contain a sizable exotic admixture of 
a tetraquark (qqqq) state. 



In an extension of the work of ( Melnitchouk et al. 



2003 Nemoto et al. 2003) in the quenched approxi 



mation, ([Takahashi and Oka 2010) have studied the 



A* (1405) on Nf = 2 CP-PACS lattices and essentially 



reach the same conclusion as ( Melnitchouk et al. 2003 
Nemoto et aL[ 12003) that the A* (1405) can not be repro- 



duced using standard baryon octet and singlet interpo- 



lating operators. 



C. Results with dynamical light and strange quarks 

The first large scale computation of the light hadron 
spectrum with a pair of light and one strange sea quark 



was performed by the MILC collaboration (Aubin et al. 
20041 iBernard et al.\ 12001 P^ With asqtad fermions on a 



one-loop Symanzik improved gauge action, they reached 
Goldstone (i.e. taste pseudoscalar) pion masses down to 
A/tt ^ 260 MeV on lattices of spatial size L ^ 2.4 fm 
and L ^ 3.4 fm at two values of the lattice spacing 
a ~ 0.09 fm and a ~ 0.12 fm. Finite volume effects were 
explicitly checked for and found to be under control. The 
fermion update algorithm used was the R-algorithm and 
explicit checks for the absence of step size dependent ef- 
fects were performed. The scale was set via b-meson 
spectroscopy, in particular the T IP-IS mass splitting, 
and physical light and strange quark masses were de- 
fined by Afjr and Mk. Ground state meson and some 
baryon masses were computed as well as the radially ex- 
cited pseudoscalar meson state. The extrapolation to 
physical pion masses was performed using various heavy 
baryon xPT motivated fit functions and a continuum ex- 
trapolation was done using g^a^ terms were possible. An 
update of these results including data from finer lattices 



as well as a comprehensive review is available in (Baza- 



vov et al. 2010a I. The resulting light hadron spectrum 



is displayed in fig. [19] Note that due to the particular 
difficulties in extracting baryon masses in the staggered 



formulation (cf. sect. Ill ) there are only predictions for a 



subset of the ground state baryons. Again, the numbers 
turn out to be in good agreement with experiment. 

A subset of the MILC ensembles with a ^ 0.12 fm and 
with a smallest pion mass of ~ 290 MeV has been studied 



J Q in (Walker-Loud et al. . 2009 ) in a mixed action setup with 



domain wall valence quarks. Comparing different chiral 
fit forms for the nucleon mass it was demonstrated that 
a simple linear fit in M^^ gives a good description of the 
given data set and extrapolates to the correct value at 
the physical point. In the same paper, this feature has 
also been found in other collaborations data. 

The PACS-CS collaboration has published results for 
the light hadron spectrum using both a chiral extrapo- 



the physical point (|Aoki et al 
2 



lation (Aoki et al. 2009a) and a direct reweighting to 

In both cases 



2010) 



N 



f 



1 nonperturbatively 0{a) improved cover 
fermions on an Iwasaki gauge action were used at a sin- 
gle lattice spacing a ~ 0.09 fm and a spatial lattice ex- 
tent of i ^ 2.9 fm. Pion masses down to ^ 150 MeV 



See also (! Davies et al.\ |2004[ l where the effects of unquenching 
are discussed for observables beyond the light hadron spectrum. 
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FIG. 19 Comparison of the Nf = 2 + 1 light hadron spectrum 



results from the MILC collaboration (Bazavov et al. 2010a I 
with experiment. Diamonds are input quantities while cir- 
cles are predictions. Experimental masses of hadrons from 
( [Amsler et aL[[2008[ ) are indicated by squares. Note that the 
plot also includes charmonium and bottomonium masses with 
some of the later ones used to set the scale. Plot reproduced 
with friendly permission of the MILC collaboration. 



were directly simulated and a reweighting to the phys- 
ical point was carried out with the lightest ensemble. 
In the extrapolated ensemble finite size effects on the 
pseudoscalar masses were corrected using SU(2) xPT at 
NLO. The tiny chiral extrapolation was performed lin- 
early in the light quark mass and Mq was used to set 
the scale. More involved chiral forms were subsequently 
investigated in (Ishikawa et al.\ 2009). Similarly in the 



reweighted ensemble the masses of the tt, K and J7 were 
used to tune to the physical point. The final result from 
the extrapolation method is plotted in fig. [20] Very simi- 
lar results have been found with the reweighting method 



as detailed in (Aoki et al. 2010). 



Full control over all systematic uncertainties at the few 
percent level was achieved in the light hadron spectrum 
calculation of the Budapest-Marseille- Wuppertal collabo- 
ration ( Durr et al. 2008 1 . They used tree level improved 
6-step stout smeared Nf = 2 -|- 1 clover fermions on a 
tree level Symanzik improved gauge action on lattices of 
spatial extent of L ~ 2.0 — 4.1 fm. Both the gauge and 
the fermion action are known to be in the correct uni- 
versality classes and the updating algorithm is exact and 
free of possible ergodicity problems. Pion masses down 
to 190 McV and three lattice spacings a ^ 0.065 fm, 
a ^ 0.85 fm and a ~ 0.125 fm were used which allowed 
for a fully controlled extrapolation to the continuum and 
the physical point with various ansatze for both. Possi- 
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I ' 1 i 1 2 

lattice hadron masses of ( Durr et al. , 2008 ) at physical M^ — 



M^/2 in physical units. The scale setting variable Mn and 
the nucleon mass are plotted vs. the square of the pion mass 
together with a fit of the data at every lattice spacing. The 
vertical dashed line represents the physical pion mass. 



ble contamination of the propagators from excited states 
were accounted for by varying the fit range. Finite vol- 
ume corrections were applied including energy shifts for 



resonant states (as described in sect. IV. C. 2 1 that allowed 
for a detailed treatment of resonant states, too. The con- 
tinuum extrapolation was performed with a term linear 
in a or a^ and chiral fits were done with both Taylor 
and NLO heavy baryon xPT with a free coefficient (see 
fig. [21] for an example extrapolation to the physical point 
and continuum limit). The above procedure allowed for a 
fully controlled calculation of the systematic uncertainty 
via the spread of the results of all analyses weighted by 
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FIG. 23 Chiral behavior of the ratio of individ- 

ual octet masses over the average octet mass Xn = 
I (M]v + Me + Afs) vs. the ratio of the square of the pion 
mass over the square average of pseudoscalar meson masses 
X^ = I (2A'/i + M^) as obtained by the QCDSF-UKQCD 



the fit quality. The ground state light hadron spectrum collabo ration. The plot is reproduced from ( [Bietenholz et al. 



was reproduced at the percent level (cf. fig. 22 ) 



2010a I with friendly permission of the QCDSF-UKQCD col- 



The QCDSF-UKQCD collaboration has recently pro- 
posed a different approach to the physical point start- 
ing from an SU(3) symmetric theory and systematicaUy 
expanding in the SU(3) breaking parameter while keep- 
M^ constant (Bietenholz et al. 2010a|b 



mg 



2M|. 



2011 ). Preliminary results at a single lattice spacing a ■ 
0.076 fm and a spatial lattice extent of L ^ 1.2 — 2.5 fm 
are displayed in fig. |23[ They show a linear depen- 
dence of the octet and decuplet masses considered and a 
good agreement with the experimentally observed hadron 
spectrum. An Nf =2 + 1 nonperturbatively improved 
single step stout smeared clover action on a tree level 
Symanzik improved gauge action was used for this study. 
Finite size corrections are not yet applied. 

There is also an ongoing effort to compute ground 
state baryons with twisted mass fermions including a 
dynamical strange quark. As the twisted mass formal- 
ism necessitates an even number of fermion flavors (cf. 
sect. II.D.3I, these 
quark {Nf = 2 



laboration. 



were found. 

The Hadron Spectrum Collaboration is using 
anisotropic lattices in order to obtain a fine time reso- 
lution of the propagators. These ensembles are mainly 
used to extract the highly excited baryon spectrum. 
The lattice spacing in time direction is tuned to be 
smaller by a factor of ^ ~ 3.5 than the lattice spacing 
in the spatial directions (Edwards et al. 2008). In 



their excited state spectroscopy studies (Bulava et al. 



2010 Dudek et al. 2011 Lin et al. 2009) they employ 



calculations also include a charm 
1 + 1). First preliminary results of 



this effort are reported in (Drach et al. 20101. 



The RBC-UKQCD collaboration has recently per- 
formed a pioneering calculation of the 77 and ry' masses 
using Nf = 2 + 1 flavor domain wall ensembles on an 
Iwasaki gauge action (Christ et al. 2010). Three pion 



Nf = 2 + 1 anisotropic clover fermions on a tree level 
tadpole improved Symanzik gauge action. A single 
spatial lattice spacing Os ~ 0.12 fm and three pion 
masses in the range 390 — 530 MeV are used. The 
scale is set with Mq. A variational method based on a 
large number (6-10) of specifically tailored interpolating 
operators are used to extract the tower of excited states 
in the different channels. Results are reported at three 
different pion masses and show a nice overall qualitative 
agreement with the experimentally observed excited 
hadron spectrum (see fig. 24). The authors emphasize 



masses in the range 400 — 700 MeV with a single lat- 
tice spacing a ~ 0.11 fm on latices with a spatial extent 
of L ^ 1.8 fm were used. A two operator basis with 
gauge fixed wall sources was used to extract the corre- 
lation functions. A mixing angle of 8 = —9.2(4.7)° and 
masses M^ = 583(15) MeV and A^ = 853(123) MeV ([Mathur et al] [2010 ) 



the need for multi hadron interpolating operators in or- 
der to reliably identify scattering states. More recently, 
also the spins of nucleon and A excitations up to spin 



7/2 have been identified by (Edwards et al. 2011). 



Ground and excited state meson spectra are also be- 
ing studied with overlap valence on dynamical domain 
wall fermions. Some preliminary results can be found in 
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FIG. 24 Comparison of part of the excited state spectrum of nucleon (left) A (middle) and fl (right) type baryons as computed 
by the hadron spectrum collaboration at three different pion masses with experiment. More details can be found in the original 
paper. The plot is reproduced from ([Bulava et al.\ |2010[) with friendly permission of the hadron spectrum collaboration. 



The quenched studies of (Mahbub et al. 2010a|c 



2009b I on the excited baryon spectrum, especially the 



excited states of the nucleon, were recently extended 
to Nf = 2 + 1 dynamical configurations by (Mahbub 
et al.\ |2010b|d[ ). Fat link irrelevant clover (FLIC) va- 
lence fermions were used on the PACS-CS dynamical en- 
sembles discussed above. Large operator bases of up to 
8 were used and signals for up to 3 excited states were 
identified. The chiral behavior of both positive and neg- 
ative nucleon excitations was studied and some evidence 
was found for the correct ordering of the negative parity 
ground state and the Roper resonance as one approaches 
physical pion masses. 



VI. CONCLUDING REMARKS 

Although it has taken over 30 years from the formula- 
tion of QCD as the theory of the strong force and Wilsons 
lattice regularization, it is fair to say that today we have 
a firm, quantitative understanding of the most relevant 
part of its particle content. It has taken so long to reach 
this level of understanding because low energy QCD is a 
very rich and nonperturbative theory. The mechanism of 
permanent quark confinement and the subsequent emer- 
gence of a particle spectrum that does not at all reflect 
the fundamental degrees of freedom required the develop- 
ment of an entirely new set of techniques that have now 
matured to a point where the experimentally observed 
spectrum of ground state, light non-singlet hadrons can 
be reproduced to an accuracy of a few percent. 

This quantitative understanding was gained in a pro- 
cess that spanned several decades. Although the funda- 
mental theory and the general strategy towards its non- 
perturbative first-principles solution was clear from the 
beginning, it required a substantial amount of conceptual 
development and physical insight. 

It is however still not a trivial task today to obtain a 



precise prediction with fully controlled uncertainties from 
QCD in the regime where it is a strongly coupled gauge 
theory. One needs to be careful of optimizing all aspects 
of the calculation to such a degree that no single one 
of them does fully dominate the total error while at the 
same time keeping the formalism simple and transpar- 
ent enough that computations are manageable in a rea- 
sonable amount of time. While ground state non-singlet 
hadron masses can be computed to a few percent accu- 
racy today, reaching the same level of precision for ex- 
cited states or singlet hadrons is still a challenging task. 
There has been substantial progress regarding the extrac- 
tion of excited states and disconnected diagram contribu- 
tions and the current understanding is approaching the 
precision level. A detailed treatment of resonant finite 
volume effects, the continuum extrapolation and even 
reaching the physical point is work currently in progress. 



Lattice calculations of the ground state, non-singlet 
hadron masses are currently trying to enter the sub- 
percent level precision region. In order to reach this goal, 
the next challenges involve a first principles treatment of 
electromagnetic and isospin breaking effects as well as an 
improved treatment of finite volume effects in the case of 
resonant states. 



In spite of these many open questions and future chal- 
lenges, we do believe however that the percent level un- 
derstanding of relevant parts of the light hadron spec- 
trum with fully controlled systematic uncertainties that 
has been achieved by lattice QCD is a milestone that 
marks the overall maturity of the method. While a lot of 
interesting problems such as excited state spectroscopy 
still require substantial work, lattice QCD today repre- 
sents a reliable tool of extracting from first principles 
properties of a strongly coupled quantum field theory. 
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